A  DYNAMIC  FRAMEWORK  FOR  SUBDIVISION  SURFACES 


By 
CHHANDOMAY  MANDAL 


A  DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 

OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 

OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 

DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 
1998 


©Copyright  1998 

by 

Chhandomay  Mandal 


To 

My  Parents, 

and  My  Brother 


ACKNOWLEDGMENTS 

I  am  deeply  grateful  to  my  advisor  Dr.  Baba  C.  Vemuri,  whose  guidance,  en- 
couragement and  support  made  the  completion  of  this  work  possible.  I  thank  him 
for  his  academic  as  well  as  non-academic  advice  and  assistance  during  my  study  at 
the  University  of  Florida,  and  for  his  impressive  example  of  integrity  and  academic 
responsibility.  I  wish  to  express  sincere  appreciation  to  my  co-advisor,  Dr.  Hong  Qin, 
for  introducing  me  to  the  wonders  of  computer  graphics,  for  patiently  guiding  me  all 
the  way  till  the  end,  and  for  assisting  me  in  various  ways.  I  would  also  like  to  thank 
Dr.  Sartaj  Sahni,  Dr.  Paul  A.  Fishwick  and  Dr.  Douglas  A.  Cenzer  for  serving  in 
my  supervisory  committee  and  advising  me  on  various  aspects  of  this  dissertation. 

My  dissertation  research  was  supported  in  part  by  the  NSF  grant  ECS-9210648 
and  the  NIH  grant  RO1-LM05944  of  Dr.  Vemuri,  and  by  the  NSF  CAREER  award 
CCR-9702103  and  DMI-9700129  of  Dr.  Qin.  I  wish  to  acknowledge  Dr.  Tim  Mcln- 
erney,  Dr.  Gregoire  Malandain,  Dr.  Hughes  Hoppe,  Dr.  Kari  Pulli,  Dr.  Dimitry 
Goldgof,  Dr.  Christina  Leonard  and  Dr.  Shang-Hong  Lai  for  helping  me  at  various 
stages  of  the  work  by  providing  various  data  sets,  software  and  figures. 

I  would  like  to  thank  the  faculty,  staff  and  students  of  the  computer  and  in- 
formation science  and  engineering  department,  who  were  always  there  to  help  me. 


In  my  years  with  the  computational  vision,  graphics  and  medical  imaging  group,  I 
have  enjoyed  working  with  a  set  of  bright  students,  including  Yanlin  Guo,  Yi  Cao, 
Li  Chen,  Li  Wang,  Shuangying  Huang,  Chinar  Kapoor,  Fengting  Chen,  Arun  Srini- 
vasan,  Jundong  Liu  and  Jun  Ye.  Special  thanks  goes  to  my  long  time  roommates, 
Raja  Chatterjee  and  Kingshuk  Majumdar,  for  helping  me  in  various  ways  during  my 
stay  at  Gainesville. 

Finally,  I  would  like  to  thank  my  parents  and  my  brother  for  their  constant  en- 
couragement and  support,  for  cheering  me  up  during  the  hard  days  and  for  providing 
an  ever  increasing  incentive  to  finish  my  graduate  study.  I  would  have  been  nowhere 
near  my  goals  without  them. 


TABLE  OF  CONTENTS 

ACKNOWLEDGMENTS    iv 

LIST  OF  FIGURES ix 

ABSTRACT    xiii 

CHAPTERS 

1  INTRODUCTION 1 

1.1  Problem  Statement 2 

1.2  Proposed  Solution 3 

1.3  Contributions 5 

1.4  Outline  of  Dissertation 8 

2  BACKGROUND 11 

2.1  Subdivision  Surfaces    11 

2.2  Physics-based  Deformable  Surface  Models     17 

2.3  Shape  Modeling  Using  Physics-based  Subdivision-surface  Model    ...  20 

2.4  Shape  and  Motion  Estimation  Using  Physics-based  Subdivision-surface 
Model 22 

2.5  Unified  Framework  for  Shape  Recovery  and  Shape  Modeling 24 

3  DYNAMIC  CATMULL-CLARK  SUBDIVISION  SURFACES     27 

3.1  Overview  of  the  Catmull-Clark  Subdivision  Scheme 27 

3.2  Formulation 30 

3.2.1  Assigning  Patches  to  Regular  Faces 31 

3.2.2  Assigning  Patches  to  Irregular  Faces     33 

3.2.3  Kinematics     36 

3.2.4  Dynamics 44 

3.2.5  Multilevel  Dynamics 46 

3.3  Finite  Element  Implementation 49 

3.3.1  Data  Structures 49 

3.3.2  Normal  Elements 51 

3.3.3  Special  Elements 51 


3.3.4  Force  Application 52 

3.3.5  Discrete  Dynamic  Equation 54 

3.3.6  Model  Subdivision    55 

3.4     Generalization  of  the  Approach 56 

4  DYNAMIC  BUTTERFLY  SUBDIVISION  SURFACES 57 

4.1  Overview  of  the  (Modified)  Butterfly  Subdivision  Scheme 57 

4.2  Formulation 61 

4.2.1  Local  Parameterization 62 

4.2.2  Dynamics 72 

4.2.3  Multilevel  Dynamics 74 

4.3  Finite  Element  Implementation 76 

4.3.1  Elemental  Mass  and  Damping  matrices 77 

4.3.2  Elemental  Stiffness  Matrix 79 

4.3.3  Other  Implementation  Issues 81 

4.4  Generalization  of  the  Approach 82 

5  UNIFIED  DYNAMIC  FRAMEWORK  FOR  SUBDIVISION  SURFACES  .  83 

5.1  Overview  of  the  Unified  Approach 83 

5.2  Unified  Dynamic  Framework  for  Catmull-Clark  Subdivision  Surfaces  .  85 

5.2.1  Local  Parameterization 88 

5.2.2  Dynamics 92 

5.2.3  Finite  Element  Implementation 92 

5.3  Unified  Dynamic  Framework  for  Loop  Subdivision  Surfaces 95 

5.3.1  Local  Parameterization 97 

5.3.2  Dynamics 100 

5.3.3  Finite  Element  Implementation 100 

5.4  A  General  Outline  of  the  Framework  for  Approximating  Subdivision 
Schemes 101 

5.5  A  General  Outline  of  the  Framework  for  Interpolatory  Subdivision 
Schemes 102 

6  MULTIRESOLUTION  DYNAMICS    103 

6.1  Overview  of  Multiresolution  Analysis  and  Wavelets 109 

6.2  Multiresolution  Analysis  for  Surfaces  of  Arbitrary  Topology     114 

6.2.1  Nested  Spaces  through  Subdivision 116 

6.2.2  Inner  Product 118 

6.2.3  Biorthogonal  Surface  Wavelets  on  Arbitrary  Manifold 119 

6.3  Multiresolution  Representation     125 

6.4  Dynamics 128 

6.5  Implementation  Details 129 


7  SYSTEM  ARCHITECTURE    133 

7.1  Topological  Information  Processing  Module     135 

7.2  Subdivision  Module 135 

7.3  Finite  Element  Analysis  Module 135 

7.4  Force  Synthesis  Module 136 

7.5  Update  Engine    136 

7.6  Display  Module 137 

8  APPLICATIONS 138 

8.1  Geometric  Modeling    138 

8.2  Shape  Recovery  from  Range  and  Volume  Data 142 

8.3  Non-rigid  Motion  Estimation 148 

8.4  Multiresolution  Editing 150 

8.5  Natural  Terrain  Modeling    151 

9  CONCLUSIONS  AND  FUTURE  WORK 155 

9.1  Conclusions 155 

9.2  Future  Directions 156 

9.2.1  Automatic  Change  of  Topology 156 

9.2.2  Local  Refinement 157 

9.2.3  Automatic  Model  Parameter  Selection    157 

9.2.4  Constraint  Imposition    157 

9.2.5  Recovery  of  Sharp  Features     158 

9.2.6  Automatic  Model  Initialization    158 

9.2.7  Improved  Synthesized  Force  Fields    158 

REFERENCES 159 

BIOGRAPHICAL  SKETCH 168 


LIST  OF  FIGURES 


2.1     Refinement  of  an  initial  control  mesh  to  obtain  the  limit  surface. 

(Courtesy  :  H.  Hoppe) 12 

3.1  Catmull-Clark  subdivision  :  (a)  initial  mesh,  (b)  mesh  obtained  after 
one  step  of  Catmull-Clark  subdivision,  and  (c)  mesh  obtained  after 
another  subdivision  step 28 

3.2  A  rectangular  mesh  and  its  limit  surface  consisting  of  4  bicubic  surface 
patches 32 

3.3  A  mesh  with  an  extraordinary  point  of  degree  3  and  its  limit  surface.       34 

3.4  Local  subdivision  around  the  extraordinary  point  and  the  limit  surface.   35 

3.5  Local  subdivision  around  the  extraordinary  point  and  the  correspond- 
ing patches  in  the  limit  surface  from  different  levels  of  subdivision.  .  .      43 

3.6  The  data  structure  used  for  dynamic  subdivision  surface  implementation.  50 

3.7  Model  subdivision  to  increase  the  degrees  of  freedom  :  (a)  evolution 
of  the  initial  mesh,  and  (b)  the  corresponding  limit  surface  evolution 
perceived  by  the  user 55 

4.1  (a)  The  control  polygon  with  triangular  faces;  (b)  mesh  obtained  after 

one  subdivision  step  using  butterfly  subdivision  rules 58 

4.2  (a)  The  contributing  vertices  in  the  j-th  level  for  the  vertex  in  the  j+1- 
th  level  corresponding  to  the  edge  between  v\  and  v32:  (b)  the  weighing 
factors  for  different  vertices 58 

4.3  (a)  The  weighing  factors  of  contributing  vertex  positions  for  an  edge 
connecting  two  vertices  of  degree  6;  (b)  the  corresponding  case  when 

one  vertex  is  of  degree  n  and  the  other  is  of  degree  6 60 

4.4  The  smoothing  effect  of  the  subdivision  process  on  the  triangles  of  the 
initial  mesh 62 


4.5  Tracking  a  point  x  through  various  levels  of  subdivision  :  (a)  initial 
mesh,  (b)  the  selected  section  (enclosed  by  dotted  lines)  of  the  mesh 
in  (a),  after  one  subdivision  step,  (c)  the  selected  section  of  the  mesh 

in  (b),  after  another  subdivision  step 64 

4.6  Different  set  of  control  points  at  a  subdivided  level  obtained  by  apply- 
ing different  subdivision  matrices  on  a  given  set  of  control  points  in  a 
coarser  mesh 67-68 

4.7  (a)  An  initial  mesh,  and  (b)  the  corresponding  limit  surface.  The  do- 
mains of  the  shaded  elements  in  the  limit  surface  are  the  corresponding 
triangular  faces  in  the  initial  mesh.  The  encircled  vertices  in  (a)  are 

the  degrees  of  freedom  for  the  corresponding  element 75 

5.1  A  control  mesh  with  an  extraordinary  vertex  of  degree  5  and  the  corre- 
sponding limit  surface  :  (a)  using  the  concepts  developed  in  Chapter  3, 
where  the  limit  surface  consists  of  quadrilateral  normal  elements  and 
a  pentagonal  special  element;  (b)  using  the  unified  concept  developed 
in  this  chapter,  where  the  limit  surface  consists  of  one  single  type  of 
quadrilateral  finite  element 86 

5.2  In  Catmull-Clark  subdivision,  each  non-boundary  quadrilateral  face  in 
the  control  mesh  has  a  corresponding  quadrilateral  patch  in  the  limit 
surface  :  (a)  control  mesh,  (b)  limit  surface 87 

5.3  (a)  The  marked  16  control  vertices  define  the  shaded  quadrilateral 
patch  associated  with  the  shaded  regular  face  in  the  control  mesh;  (b) 
the  marked  14  control  vertices  define  the  shaded  quadrilateral  patch 
associated  with  the  shaded  irregular  face  in  the  control  mesh 91 

5.4  (a)  The  control  polygon  with  triangular  faces;  (b)  mesh  obtained  after 

one  subdivision  step  using  Loop's  subdivision  rules 96 

5.5  An  initial  mesh  and  the  corresponding  limit  surface  obtained  using 
Loop's  subdivision  rules.  The  domains  of  the  shaded  triangular  patches 
in  the  limit  surface  are  the  corresponding  triangular  faces  in  the  initial 
mesh.  The  encircled  vertices  are  the  control  vertices  for  the  corre- 
sponding triangular  patch  in  the  limit  surface 98 

5.6  Each  triangular  patch  in  the  limit  surface  can  be  associated  with  a 
non-boundary  triangular  face  in  the  initial  mesh,  which  in  turn  can  be 
parameterized  over  a  triangle  with  vertices  at  (0,0),  (1,0)  and  (0, 1).  .      99 

6.1     Schematic  block  diagram  of  the  multilevel  dynamics  approach 104 


6.2  Multilevel  dynamics  vs.  multiresolution  dynamics 106 

6.3  Representation  of  the  degrees  of  freedom  in  multilevel  dynamics  and 
multiresolution  dynamics  approach 108 

6.4  The  filter  bank 112 

6.5  Subdivision  refinement  of  a  tetrahedron 119 

6.6  Wavelet  construction 123 

7.1     System  architecture 134 

8.1  (a),  (b),  (c)  and  (d)  :  Initial  shapes  (obtained  applying  Catmull-Clark 
subdivision  rules  on  control  meshes);  (e),  (f),  (g)  and  (h)  :  the  corre- 
sponding modified  shapes  after  interactive  force  application 139 

8.2  (a),  (b),  (c)  and  (d)  :  Initial  shapes  (obtained  applying  modified  but- 
terfly subdivision  rules  on  control  meshes);  (e),  (f),  (g)  and  (h)  :  the 
corresponding  modified  shapes  after  interactive  sculpting  via  force  ap- 
plication  141 

8.3  (a)  Range  data  of  a  bulb  along  with  the  initialized  model,  (b)  an 
intermediate  stage  of  evolution,  and  (c)  the  fitted  dynamic  Catmull- 
Clark  subdivision  surface  model 143 

8.4  (a)  Range  data  of  an  anvil  along  with  the  initialized  model,  (b)  an 
intermediate  stage  of  evolution,  and  (c)  the  fitted  dynamic  Catmull- 
Clark  subdivision  surface  model 144 

8.5  (a)  Range  data  of  a  head  along  with  the  initialized  model,  (d)  the 
fitted  dynamic  butterfly  subdivision  model,  and  (c)  visualization  of 

the  shape  from  another  view  point 144 

8.6  (a)  A  slice  from  a  brain  MRI,  (b)  initialized  model  inside  the  region 
of  interest  superimposed  on  the  slice,  (c)  the  fitted  model  superim- 
posed on  the  slice,  and  (d)  a  3D  view  of  the  dynamic  Catmull-Clark 
subdivision  surface  model  fitted  to  the  cerebellum 146 

8.7  (a)  Data  points  identifying  the  boundary  of  the  caudate  nucleus  on 
a  MRI  slice  of  human  brain,  (b)  data  points  (from  all  slices)  in  3D 
along  with  the  initialized  model,  and  (c)  the  fitted  dynamic  butterfly 
subdivision  model 146 


.8     Snapshots  from  the  animation  of  canine  heart  motion  over  a  cardiac 

cycle  using  the  dynamic  butterfly  subdivision  model 149 

.9  (a)  The  initialized  model  along  with  the  data  set;  (b)  the  fitted  model 
with  two  subdivisions  on  the  initial  mesh,  along  with  attached  springs 
for  editing.  The  model  after  editing  (c)  at  lower  resolution,  (d)  at  the 
same  resolution  of  the  fitted  model,  and  (e)  at  higher  resolution.  .  .  .    150 

1.10  (a)  Discrete  elevation  data  set  (4096  points),  (b)  fitted  dynamic  but- 
terfly subdivision  surface  model  with  841  vertices  (without  noise  ad- 
dition), and  (c)  fitted  dynamic  subdivision  surface  model  with  841 
vertices  when  Gaussian  noise  is  added 153 

i.ll  Synthesized  natural  terrain  of  different  roughness  using  the  dynamic 
butterfly  subdivision  surface  model  with  841  vertices  from  a  data  set 
of  3099  elevation  values 154 


Abstract  of  Dissertation 

Presented  to  the  Graduate  School  of  the  University  of  Florida 

in  Partial  Fulfillment  of  the  Requirements  for  the 

Degree  of  Doctor  of  Philosophy 


A  DYNAMIC  FRAMEWORK  FOR  SUBDIVISION  SURFACES 

By 

Chhandomay  Mandal 

December  1998 


Chairman:  Dr.  Baba  C.  Vemuri 

Cochairman  :  Dr.  Hong  Qin 

Major  Department:  Computer  and  Information  Science  and  Engineering 


Subdivision  surfaces  are  extensively  used  to  model  smooth  shapes  of  arbitrary 
topology.  Recursive  subdivision  on  a  user-defined  initial  control  mesh  generates  a 
visually  pleasing  smooth  surface  in  the  limit.  However,  users  have  to  carefully  select 
the  initial  mesh  and/or  manipulate  the  control  vertex  positions  at  different  levels 
of  subdivision  hierarchy  to  satisfy  the  functional  and  aesthetic  requirements  in  the 
smooth  limit  surface.  This  modeling  drawback  results  from  the  lack  of  direct  manip- 
ulation tools  for  the  limit  surface.  In  this  dissertation,  techniques  from  physics-based 
modeling  are  integrated  with  geometric  subdivision  methodology,  and  a  dynamic 
framework  is  presented  for  direct  manipulation  of  the  smooth  limit  surface  generated 
by  the  subdivision  schemes  using  physics-based  "force"  tools. 


In  a  typical  subdivision  scheme,  the  user  starts  with  an  initial  control  mesh 
which  is  refined  recursively  using  a  fixed  set  of  subdivision  rules,  and  a  smooth  sur- 
face is  produced  in  the  limit.  Most  often  this  limit  surface  does  not  have  an  analytic 
expression,  and  hence  poses  challenging  problems  in  incorporating  mass  and  damp- 
ing distribution  functions,  internal  deformation  energy,  forces,  and  other  physical 
quantities  required  to  develop  a  physics-based  subdivision  surface  model.  In  this 
dissertation,  local  parameterization  techniques  suitable  for  embedding  the  geometric 
subdivision  surface  model  in  a  physics-based  modeling  framework  have  been  devel- 
oped. Specific  local  parameterization  techniques  have  been  fully  developed  for  the 
Catmull-Clark,  modified  butterfly  and  the  Loop  subdivision  schemes.  Techniques 
for  assigning  material  properties  to  geometric  subdivision  surfaces  are  derived,  and 
a  motion  equation  for  the  dynamic  model  has  been  formulated  using  Lagrangian 
dynamics.  Furthermore,  advantages  of  the  physics-based  deformable  models  are  in- 
corporated into  the  conventional  subdivision  schemes,  and  a  dynamic  hierarchical 
control  of  this  model  is  introduced.  Finally,  a  multiresolution  representation  of  the 
control  mesh  is  developed  and  a  unified  approach  for  deriving  subdivision  surface- 
based  finite  elements  is  presented. 

The  proposed  dynamic  framework  enhances  the  applicability  of  the  subdivision 
surfaces  in  modeling  applications.  It  is  also  useful  for  hierarchical  shape  recovery 
from  large  range  and  volume  data  sets,  as  well  as  for  non-rigid  motion  tracking  from 
a  temporal  sequence  of  data  sets.  Multiresolution  representation  of  the  initial  mesh 
controlling  the  smooth  limit  surface  enables  global  and  local  editing  of  the  shape  as 


desired  by  the  modeler.  This  dynamic  framework  has  also  been  used  for  synthesizing 
natural  terrain  models  from  sparse  elevation  data. 


CHAPTER  1 
INTRODUCTION 


Generating  smooth  surfaces  of  arbitrary  topology  poses  a  grand  challenge  for 
the  computer  graphics,  computer-aided  geometric  design  and  scientific  visualization 
researchers.  The  existing  techniques  for  modeling  smooth  complex  shapes  can  be 
broadly  classified  into  two  distinct  categories  namely,  (1)  explicit  patching  using 
parametric  surfaces  and  (2)  subdivision  surfaces. 

The  explicit  patching  technique  involves  partitioning  the  model  surface  into  a 
collection  of  individual  parametric  surface  patches.  Adjacent  surface  patches  are 
then  explicitly  stitched  together  using  continuity  constraints.  This  explicit  stitching 
process  is  very  complicated  for  modeling  smooth  surfaces  of  arbitrary  topology  due 
to  the  continuity  constraints  which  need  to  be  satisfied  along  the  patch  boundaries. 

Subdivision  surfaces  are  simple  procedural  models  which  offer  an  alternate  rep- 
resentation for  the  smooth  surfaces  of  arbitrary  topology.  A  typical  recursive  subdi- 
vision scheme  produces  a  visually  pleasing  smooth  surface  in  the  limit  by  repeated 
application  of  a  fixed  set  of  refinement  rules  on  an  user-defined  initial  control  mesh. 
The  initial  control  mesh  is  a  simple  polygonal  mesh  of  the  same  topological  type 
as  that  of  the  smooth  surface  to  be  modeled.  At  each  step  of  subdivision,  a  finer 
polygonal  mesh  with  more  vertices  and  faces  is  constructed  from  the  previous  one  via 


the  refinement  process,  and  the  smooth  surface  is  obtained  in  the  limit.  However,  a 
few  subdivision  steps  on  the  initial  mesh  generally  suffice  to  approximate  the  smooth 
surface  for  all  practical  purposes.  Various  sets  of  rules  lead  to  different  subdivision 
schemes  which  mainly  differ  in  the  smoothness  property  of  the  resulting  limit  surface 
and/or  the  type  of  initial  mesh  (i.e.  triangular,  quadrilateral  etc.)  chosen. 

1.1     Problem  Statement 

To  model  a  smooth  surface  of  arbitrary  topology  using  subdivision  surfaces,  first 
a  polygonal  mesh  of  same  topology  needs  to  be  chosen  as  the  initial  mesh.  This  ini- 
tial mesh,  also  known  as  the  control  mesh,  is  refined  via  recursive  subdivision  using 
a  fixed  set  of  rules,  and  a  smooth  surface  of  the  desired  topology  is  obtained  in  the 
limit.  However,  users  have  to  carefully  select  the  initial  mesh  and/or  manipulate 
the  control  vertex  positions  at  different  levels  of  subdivision  hierarchy  in  order  to 
satisfy  the  functional  and  aesthetic  requirements  on  the  smooth  limit  surface.  For 
example,  to  obtain  a  desired  effect  on  the  smooth  limit  surface,  it  might  be  necessary 
to  reposition  a  handful  of  vertices  in  the  mesh  obtained  after  one  subdivision  step, 
or  a  large  number  of  vertices  might  need  to  be  moved  in  the  mesh  produced  after 
three  subdivision  steps!  This  process  is  not  intuitive  and  at  best  laborious.  Despite 
the  presence  of  a  variety  of  subdivision  schemes  in  the  computer  graphics  and  geo- 
metric modeling  literature,  there  is  no  direct  and  natural  way  of  manipulating  the 
limit  surface.  The  current  state-of-the-art  only  permits  the  modeler  to  interactively 
obtain  the  desired  effects  on  the  smooth  surface  by  kinematically  manipulating  the 
vertex  positions  at  various  levels  of  the  subdivision  hierarchy.  In  this  dissertation,  tlie 


challenging  problem  of  directly  manipulating  the  smooth  limit  surface  at  arbitrary 
locations/areas  is  addressed  and  a  novel  solution  is  presented  where  the  modeler  can 
not  only  manipulate  the  smooth  limit  surface  directly  but  also  control  the  extent  of 
the  manipulation  effect  (i.e.,  global  or  local  manipulation)  on  the  limit  surface. 

Subdivision  surfaces  have  also  been  used  for  recovering  shapes  from  a  given  set  of 
points  in  3D.  However,  most  of  the  existing  subdivision  surface-based  shape  recovery 
techniques  resort  to  complex  algorithms  to  derive  a  mesh  for  the  underlying  shape, 
and  then  mesh  optimization  techniques  are  used  to  obtain  a  compact  representation 
of  the  same.  Nevertheless,  this  process  yields  a  control  mesh  of  the  smooth  subdi- 
vision surface  representing  the  underlying  shape  that  typically  uses  a  large  number 
of  degrees  of  freedom  (control  vertices)  for  representation.  In  this  dissertation,  an 
efficient  hierarchical  method  of  recovering  shapes  from  range  and  volume  data  sets  is 
proposed  where  the  control  mesh  of  the  smooth  limit  surface  will  use  very  few  degrees 
of  freedom  for  representation. 

1.2     Proposed  Solution 

Physics-based  modeling  techniques  offer  a  potential  solution  to  the  problem 
of  directly  manipulating  the  smooth  limit  surface  generated  by  the  recursive  sub- 
division procedure.  In  the  physics-based  modeling  paradigm,  a  deformable  model 
is  derived  by  assigning  mass,  damping,  stiffness  and  other  physical  properties  to  a 
surface  model.  The  model  is  deformed  by  applying  synthesized  forces  and  this  de- 
formation is  governed  by  physical  laws.  Now,  if  the  purely  geometric  subdivision 
surface  models  can  be  embedded  in  a  physics-based  modeling  paradigm,  then  the 


smooth  limit  surface  can  be  deformed  directly  by  using  physics-based  "force"  tools. 
However,  this  procedure-based  surface  model  obtained  through  the  subdivision  pro- 
cess does  not  have  a  closed-form  analytic  formulation  in  general,  and  hence  poses 
challenging  problems  to  incorporate  mass  and  damping  distribution  functions,  inter- 
nal deformation  energy,  forces,  and  other  physical  quantities  required  to  develop  a 
physics-based  model.  Techniques  of  locally  parameterizing  the  smooth  limit  surface 
generated  by  various  subdivision  schemes  are  proposed  in  this  dissertation.  Once  a 
suitable  local  parameterization  scheme  is  developed  for  a  specific  subdivision  scheme, 
a  dynamic  framework  is  provided  for  the  corresponding  subdivision  scheme  where 
the  modeler  can  directly  manipulate  the  smooth  limit  surface  by  using  synthesized 
forces.  At  the  same  time,  a  dynamic  framework  for  the  wavelets  derived  using  subdi- 
vision schemes  will  assist  in  adopting  a  multiresolution  representation  of  the  control 
mesh  defining  the  smooth  limit  surface,  and  the  modeler  can  control  the  extent  of 
manipulation  effect  by  choosing  a  desired  level  of  editing.  Synthesized  force  applica- 
tion at  a  lower  resolution  will  yield  a  global  effect  whereas  manipulations  at  a  finer 
resolution  will  have  localized  effects  on  the  limit  surface.  The  motion  of  this  physics- 
based  deformable  subdivision  surface  model  is  governed  by  a  second-order  differential 
equation,  which  is  solved  numerically  using  the  finite  element  method.  New  types  of 
finite  elements  for  the  chosen  subdivision  scheme  are  also  presented  for  representing 
the  smooth  limit  surface. 

The  proposed  dynamic  framework  for  the  subdivision  surfaces  provides  an  effi- 
cient solution  to  the  shape  recovery  problem  as  well.   A  simple  subdivision  surface 


mode!  with  very  few  vertices  in  the  control  mesh  can  be  initialized  (positioned)  fully 
inside  a  given  set  of  points  in  3D.  The  initialized  model  will  be  deformed  by  applying 
forces  synthesized  from  the  given  data  points.  When  an  equilibrium  is  obtained,  the 
number  of  vertices  in  the  control  mesh  can  be  increased  via  a  subdivision  step  on  the 
current  control  mesh  thereby  increasing  the  degrees  of  freedom  for  model  representa- 
tion, and  a  new  equilibrium  with  a  better  fit  to  the  given  data  set  can  be  obtained. 
This  process  can  be  repeated  till  a  prescribed  error  in  fit  is  achieved.  Similar  approach 
can  be  taken  for  shape  recovery  from  volume  data  sets  as  well  where  a  different  type 
of  synthesized  force  needs  to  be  specified.  The  hierarchical  shape  recovery  process 
ensures  a  compact  representation  of  the  recovered  smooth  limit  surface  using  very 
few  degrees  of  freedom. 

L2 Contributions 

In  this  dissertation,  techniques  from  physics-based  modeling  are  integrated  with 
geometric  subdivision  methodology  to  present  a  scheme  for  directly  manipulating 
the  smooth  limit  surface  generated  by  the  subdivision  process.  As  a  result,  unlike 
the  existing  geometric  solutions  that  only  allow  the  operations  on  control  vertices, 
the  proposed  methodology  and  algorithms  permit  the  user  to  physically  modify  the 
shape  of  subdivision  surfaces  at  desired  locations  via  application  of  forces.  This  gives 
the  user  a  "virtual"  clay /play-dough  modeling  environment.  The  proposed  model 
can  be  edited  directly  in  a  hierarchical  fashion  using  synthesized  forces.  Also,  this 
physics-based  subdivision  surface  model  efficiently  recovers  shapes  as  well  as  non- 
rigid  motions  from  large  range  and  volumetric  data  sets.  Note  that  this  dissertation 


neither  proposes  a  new  subdivision  technique  nor  provides  a  different  interpretation 
of  any  existing  subdivision  technique,  but  integrates  the  advantages  of  subdivision 
surface-based  and  physics-based  modeling  techniques  to  solve  important  theoretical 
and  practical  problems.  Although  the  principles  of  physics-based  modeling  are  well 
understood  by  the  graphics  experts  and  modeling  researchers,  this  dissertation  will 
greatly  advance  the  state  of  the  art  in  physics-based  shape  modeling  due  to  the 
contributions  listed  below. 

•  Local  parameterization  techniques  for  the  smooth  limit  surface  generated  by 
various  subdivision  schemes  are  systematically  derived  in  a  hierarchical  frame- 
work, and  subsequently  the  initial  control  polyhedron  can  be  viewed  as  the 
parametric  domain  of  the  physics-based  smooth  limit  surface. 

•  The  smooth  dynamic  subdivision  surface  in  the  limit  is  treated  as  a  "real" 
physical  object  represented  by  a  set  of  novel  finite  elements.  The  basis  (shape) 
functions  of  these  new  variety  of  finite  elements  are  derived  using  the  subdivision 
schemes.  The  proposed  finite  element  methods  will  prove  to  be  useful  not  only 
in  the  areas  of  computer  graphics  and  geometric  design,  but  also  in  engineering 
analysis  as  well. 

•  The  subdivision  techniques  are  used  to  create  a  surface  model  that  incorporates 
mass  and  damping  distribution  functions,  internal  deformation  energy,  forces, 
and  other  physical  quantities.  The  motion  equations  are  also  systematically 
derived  for  this  dynamic  subdivision  surface  model  whose  degrees  of  freedom  are 
the  collection  of  initial  user-specified  control  vertices.  Therefore,  the  advantages 


of  both  the  physics-based  modeling  philosophy  and  the  geometric  subdivision 
schemes  are  incorporated  within  a  single  unified  framework. 

•  Users  will  be  able  to  manipulate  this  physics-based  model  in  an  arbitrary  region, 
and  the  model  will  respond  naturally  (just  like  a  real-world  object  would)  to 
this  force  application.  This  shape  deformation  is  quantitatively  characterized 
by  the  time-varying  displacement  of  control  points  that  uniquely  define  the 
geometry  of  the  limit  surface. 

•  The  dynamic  framework  for  wavelets  derived  using  subdivision  schemes  en- 
ables a  multiresolution  representation  of  the  control  mesh  defining  the  smooth 
limit  surface.  This  provides  additional  flexibility  to  the  physics-based  modeling 
framework  since  the  modeler  is  free  to  choose  the  desired  editing  level.  If  a 
lower  resolution  editing  level  is  chosen,  the  synthesized  force  application  on  the 
smooth  limit  surface  will  have  global  effects,  whereas  editing  with  force-based 
tools  at  a  finer  resolution  will  yield  a  localized  effect. 

•  Algorithms  and  procedures  are  developed  which  approximate  the  proposed  new 
finite  elements  using  a  collection  of  linear  and/or  bilinear  finite  elements  subject 
to  the  implicit  geometric  constraints  enforced  by  the  subdivision  rules.  This 
hierarchically-structured  approximation  can  satisfy  any  user-specified  error  tol- 
erance. 

The  proposed  dynamic  framework  enhances  the  applicability  of  subdivision  sur- 
face models  in  various  application  areas.    It  provides  a  direct  and  intuitive  way  of 
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manipulating  shapes  in  a  hierarchical  fashion  for  geometric  modeling  applications. 
It  has  also  been  successfully  used  for  efficient  and  hierarchical  shape  recovery  from 
range  and  volume  data  sets  as  well  as  for  tracking  a  shape  of  interest  from  a  time 
sequence  of  range  or  volume  data  sets.  Finally,  the  dynamic  framework  of  subdivision 
surfaces  is  combined  with  a  variant  of  a  fractal  surface  synthesis  technique  to  present 
a  novel  natural  terrain  modeling  method. 

1.4    Outline  of  Dissertation 

The  rest  of  the  dissertation  is  organized  as  follows  : 

Chapter  2  contains  an  overview  of  the  subdivision  surfaces  along  with  a  review  of 
the  related  literature.  The  motivation  for  embedding  the  subdivision  surface  models 
in  a  physics-based  modeling  framework  is  discussed,  and  the  proposed  model  is  com- 
pared with  the  existing  physics-based  models.  The  advantages  of  a  unified  framework 
for  shape  modeling  and  shape  recovery  are  also  pointed  out  in  this  chapter. 

Chapter  3  provides  a  dynamic  framework  for  the  Catmull-Clark  subdivision 
scheme,  which  is  one  of  the  most  popular  subdivision  techniques  for  modeling  compli- 
cated objects  of  arbitrary  genus.  Analytic  formulation  of  the  limit  surface  generated 
by  the  Catmull-Clark  subdivision  scheme  is  derived  and  the  "physical"  quantities 
required  to  develop  the  dynamic  model  are  introduced  [57,  60,  76].  The  governing 
dynamic  differential  equation  is  derived  using  Lagrangian  mechanics  and  is  imple- 
mented using  a  finite  element  technique. 
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In  Chapter  4,  a  dynamic  framework  is  provided  for  the  butterfly  subdivision 
scheme,  another  popular  subdivision  technique  for  modeling  smooth  surfaces  of  arbi- 
trary topology.  A  local  parameterization  scheme  for  the  butterfly  subdivision  surface 
is  derived  in  a  hierarchical  style.  The  physics-based  butterfly  subdivision  model  is 
formulated  as  a  set  of  novel  finite  elements  which  are  optimally  approximated  by 
a  collection  of  standard  finite  elements  subjected  to  implicit  geometric  constraints 
[58,  61]. 

Chapter  5  presents  an  unified  approach  for  providing  a  dynamic  framework  for 
subdivision  surfaces  in  general.  In  particular,  it  has  been  shown  that  the  limit  surface 
obtained  using  any  subdivision  scheme  can  be  viewed  as  a  collection  of  either  quadri- 
lateral or  triangular  finite  elements  whose  basis  (shape)  functions  can  be  derived 
using  the  chosen  subdivision  scheme  [59]. 

In  Chapter  6,  a  dynamic  framework  is  presented  for  the  subdivision  surface- 
based  wavelet  schemes.  This  dynamic  framework  enables  a  multiresolution  represen- 
tation of  the  control  mesh  defining  the  smooth  limit  surface  and  consequently,  the 
modeler  can  control  the  extent  of  the  effect  of  manipulation  on  the  limit  surface  by 
choosing  a  proper  editing  level. 

In  Chapter  7,  a  system  that  integrates  the  implementation  of  the  concepts  pro- 
posed in  earlier  chapters  is  presented.  Several  modules  comprise  the  entire  system, 
and  the  functionality  of  these  modules  are  discussed  in  this  chapter. 
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The  proposed  modeling  framework  is  used  in  various  application  areas  and  the 
results  are  presented  in  Chapter  8.  The  proposed  dynamic  model  has  been  success- 
fully used  in  geometric  modeling,  shape  recovery  and  non-rigid  motion  estimation 
applications.  A  novel  technique  for  natural  terrain  modeling  is  also  presented  where 
the  dynamic  framework  is  combined  with  a  variant  of  a  fractal  surface  synthesis 
technique. 

Finally,  conclusions  are  drawn  in  Chapter  9  where  future  directions  of  research 
are  also  pointed  out. 


CHAPTER  2 

BACKGROUND 


In  this  chapter,  a  brief  overview  of  the  subdivision  surfaces  is  presented  followed 
by  a  review  of  the  previous  work  done  in  the  area  of  subdivision  surfaces.  This  is 
followed  by  an  overview  of  the  physics-based  models  along  with  a  review  of  the  re- 
lated literature.  The  motivating  factors  for  embedding  geometric  subdivision  surface 
models  in  a  physics-based  modeling  paradigm  are  presented  in  the  next  section.  The 
advantages  of  the  proposed  dynamic  models  over  the  existing  physics-based  models 
for  shape  recovery  are  discussed  and  finally  advantages  of  a  unified  framework  for 
shape  modeling  and  shape  recovery  are  also  presented. 

2J Subdivision  Surfaces 

The  input  to  any  subdivision  scheme  is  an  initial  mesh  (also  known  as  control 
mesh),  M°  =  (V°,  F°),  which  is  a  collection  of  vertices  V°  and  a  collection  of  faces 
F°.  The  subdivision  surface  S(M°)  associated  with  the  initial  mesh  M°  is  defined 
as  the  limit  of  the  recursive  application  of  the  refinement  R  as  shown  in  Fig.2.1. 
The  refinement  R,  when  applied  on  a  mesh  Mk  =  (Vk,Fk),  creates  a  refined  mesh 
M*+1  =  (V*+l,  F*+1)  where  the  vertices  in  Vk+l  are  computed  as  affine  combinations 
of  the  vertices  in  V*  and  the  faces  in  Fh+1  are  obtained  by  splitting  each  face  in  Fk 
into  a  fixed  number  of  sub-faces. 
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(a)  Contra]  mesh  M 


(b)  Mesh  M1  =  R{M) 


(c)Mesh.fcF=fi(ic(M)) 


(d)  Limit  surface  5(M)  =  Roa(M) 


Figure  2.1.  Refinement  of  an  initial  control  mesh  to  obtain  the  limit  surface.  (Cour- 
tesy :  H.  Hoppe) 
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It  may  be  noted  that  the  vertices  introduced  through  the  subdivision  process 
have  a  fixed  degree  in  genera]  (4  in  case  of  quadrilateral  meshes  and  6  in  case  of  tri- 
angular meshes).  The  vertices  which  do  not  have  this  degree  are  called  extraordinary 
vertices.  The  number  of  extraordinary  vertices  are  very  few  as  these  vertices  are  not 
introduced  via  subdivision  refinement  in  general.  For  example,  in  the  Catmull-Clark 
subdivision  scheme  (defined  on  quadrilateral  meshes)  the  number  of  extraordinary 
vertices  does  not  change  after  the  first  subdivision  on  the  initial  mesh,  whereas  in  the 
modified  butterfly  subdivision  scheme  (defined  on  triangular  meshes)  the  subdivision 
process  never  introduces  an  extraordinary  vertex.  The  limit  surface  defined  by  the 
subdivision  process  is  most  often  C1  (first  derivative  continuous)  or  C2  (second  deriva- 
tive continuous)  depending  on  the  subdivision  rules  expect  at  very  few  extraordinary 
points  corresponding  to  the  extraordinary  vertices  in  the  mesh.  Specific  subdivision 
schemes  along  with  the  properties  of  their  limit  surfaces  will  be  discussed  in  later 
chapters. 

A  lot  of  research  have  been  carried  out  on  subdivision  surfaces,  which  are  mainly 
focused  either  on  the  development  of  a  new  subdivision  technique  or  on  analyzing 
the  smoothness  properties  of  the  limit  surface  generated  by  a  subdivision  scheme. 
However,  in  this  dissertation  the  focus  is  entirely  different,  namely,  embedding  the 
subdivision  surfaces  in  a  physics-based  modeling  paradigm  to  enhance  the  applicabil- 
ity of  these  surface  models.  In  the  rest  of  this  section,  these  "traditional"  techniques 
of  subdivision  surfaces  are  reviewed. 
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Chaikin  [14]  first  introduced  the  concept  of  subdivision  to  the  computer  graphics 
community  for  generating  a  smooth  curve  from  a  given  control  polygon.  During  the 
last  two  decades,  a  wide  variety  of  subdivision  schemes  for  modeling  smooth  surfaces 
of  arbitrary  topology  have  been  derived  following  Chaikin's  pioneering  work  on  curve 
generation.  In  general,  these  subdivision  schemes  can  be  categorized  into  two  distinct 
classes  namely,  (1)  approximating  subdivision  techniques  and  (2)  interpolator  sub- 
division techniques.  The  limit  surface  obtained  using  an  approximating  subdivision 
technique  only  approximates  the  initial  mesh,  whereas  the  limit  surface  interpolates 
the  initial  mesh  in  case  of  interpolatory  subdivision  techniques. 

Among  the  approximating  schemes,  the  techniques  of  Doo  and  Sabin  [27,  83] 
generalize  the  idea  of  obtaining  uniform  biquadratic,  patches  from  a  rectangular  con- 
trol mesh.  This  scheme  is  an  example  of  vertex  subdivision  scheme  where  a  vertex 
surrounded  by  k  faces  is  split  into  k  sub-vertices,  one  for  each  face.  Doo  and  Sabin's 
subdivision  technique  can  be  applied  to  any  mesh  of  arbitrary  topology,  and  the  re- 
sulting smooth  limit  surface  would  be  biquadratic  B-splines,  except  at  extraordinary 
points  corresponding  to  the  extraordinary  vertices  in  the  control  mesh. 

Catmull  and  Clark  [10]  developed  a  method  for  recursively  generating  a  smooth 
surface  from  a  polyhedral  mesh  of  arbitrary  topology.  The  Catmull-Clark  subdivision 
surface,  defined  by  an  arbitrary  initial  mesh,  can  be  reduced  to  a  set  of  standard 
bicubic  B-spline  patches  except  at  a  finite  number  of  degenerate  points.  This  scheme 
is  an  example  of  face  subdivision  scheme,  where  a  fc-sided  face  is  split  into  k  sub- 
faces.  The  details  of  this  subdivision  scheme  are  discussed  in  Section  3.1.  Halstead 
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et  al.  [38]  proposed  an  algorithm  to  construct  a  Catmull-Clark  subdivision  surface 
that  interpolates  the  vertices  of  a  mesh  of  arbitrary  topology. 

Loop  [51]  presented  a  subdivision  scheme  based  on  the  generalization  of  quartic 
triangular  B-splines  for  triangular  meshes.  Loop's  subdivision  rules  are  presented  in 
Section  5.3.  Hoppe  et  al.  [40]  extended  Loop's  work  to  produce  piecewise  smooth 
surfaces  with  discontinuities  in  selected  areas  of  the  limit  surface.  They  introduced 
new  subdivision  rules  that  allow  for  sharp  features  such  as  creases  and  corners  in  the 
limit  surface. 

Peters  and  Reif  [73]  proposed  a  simple  subdivision  scheme  for  smoothing  poly- 
hedra.  Their  refinement  rules  yield  a  C1  surface  that  has  a  piecewise  quadratic  pa- 
rameterization except  at  isolated  extraordinary  points.  Most  recently,  non-uniform 
Doo-Sabin  and  Catmull-Clark  surfaces  that  generalize  non-uniform  tensor  product 
B-spline  surfaces  to  arbitrary  topologies  were  introduced  by  Sederberg  et  al.  [88].  All 
the  schemes  mentioned  above  generalize  recursive  subdivision  schemes  for  generating 
limit  surfaces  with  a  known  parameterization.  Various  issues  involved  with  character 
animation  using  these  approximating  subdivision  schemes  were  discussed  at  length 
by  DeRose  et  al.  [25]. 

The  most  well-known  interpolation-based  subdivision  scheme  is  the  "butterfly" 
algorithm  proposed  by  Dyn  et  al.  [30].  Butterfly  subdivision  method,  like  other  sub- 
division schemes,  makes  use  of  a  small  number  of  neighboring  vertices  for  subdivision. 
It  requires  simple  data  structures  and  is  extremely  easy  to  implement.  However,  it 
needs  a  topologically  regular  setting  of  the  initial  (control)  mesh  in  order  to  obtain 
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a  smooth  C"  limit  surface.  A  variant  of  this  scheme  with  better  smoothness  prop- 
erties can  be  found  in  Dyn  et  al.  [29].  Zorin  et  al.  [109]  developed  an  improved 
interpolatory  subdivision  scheme  (which  we  call  the  modified  butterfly  scheme)  that 
retains  the  simplicity  of  the  butterfly  scheme  and  results  in  much  smoother  surfaces 
even  from  irregular  initial  meshes.  These  interpolatory  subdivision  schemes  have 
extensive  applications  in  wavelets  on  manifolds,  multiresolution  decomposition  of 
polyhedral  surfaces  and  multiresolution  editing. 

A  variational  approach  for  interpolatory  refinement  has  been  proposed  by  Kobbelt 
[43,  44]  and  by  Kobbelt  and  Schroder  [46].  In  this  approach,  the  vertex  positions  in 
the  refined  mesh  at  each  subdivision  step  are  obtained  by  solving  an  optimization 
problem.  Therefore,  these  schemes  are  global,  i.e.,  every  new  vertex  position  depends 
on  all  the  vertex  positions  of  the  coarser  level  mesh.  The  local  refinement  property 
which  makes  the  subdivision  schemes  attractive  for  implementation  in  the  computer 
graphics  applications  is  not  retained  in  the  variational  approach. 

The  derivation  of  various  mathematical  properties  of  the  smooth  limit  surface 
generated  by  the  subdivision  algorithms  is  rather  complex.  Doo  and  Sabin  [28]  first 
analyzed  the  smoothness  behavior  of  the  limit  surface  using  the  Fourier  transform  and 
an  eigen-analysis  of  the  subdivision  matrix.  Ball  and  Storry  [3,  4]  and  Reif  [80]  further 
extended  Doo  and  Sabin's  prior  work  on  continuity  properties  of  subdivision  surfaces 
by  deriving  various  necessary  and  sufficient  conditions  on  smoothness  for  different 
subdivision  schemes.  Specific  subdivision  schemes  were  analyzed  by  Schweitzer  [87], 
Habib  and  Warren  [37],  Peters  and  Reif  [74]  and  Zorin  [108].   Most  recently,  Stain 
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[90]  presented  a  method  for  exact  evaluation  of  Catmull-Clark  subdivision  surfaces 
at  arbitrary  parameter  values.  It  may  be  noted  that  the  focus  of  this  dissertation 
is  not  on  deriving  a  new  subdivision  technique  or  analyzing  a  subdivision  technique, 
but  on  deriving  a  deformable  surface  model  using  these  subdivision  schemes. 

2.2     Phvsics-based  Deformable  Surface  Models 

Shape  models  can  be  broadly  categorized  into  two  types  -  lumped  parameter 
models  and  distributed  parameter  models.  A  lumped  parameter  model  uses  small 
number  of  parameters  to  represent  model  geometry,  whereas  a  distributed  parameter 
model  uses  large  number  of  degrees  of  freedom  for  model  representation.  An  example 
of  a  lumped  parameter  surface  model  is  a  superquadric.  Spline  surfaces  serve  as  a 
good  example  for  distributed  parameter  model. 

Deformable  surfaces  are  distributed  parameter  models  which  use  large  number 
of  degrees  of  freedom  for  representing  the  model  geometry.  A  large  variety  of  shapes 
can  be  modeled  using  this  type  of  model,  but  handling  a  large  number  of  degrees  of 
freedom  can  be  cumbersome.  However,  the  large  degrees  of  freedom  of  a  deformable 
model  are  embedded  in  a  physics-based  framework  to  allow  only  a  "physically  mean- 
ingful" model  behavior.  Various  types  of  energies  are  assigned  to  the  model  using  the 
degrees  of  freedom,  and  the  model  "deforms"  to  find  an  equilibrium  state  with  min- 
imal energy.  The  motion  equation  is  derived  using  Lagrangian  dynamics  [36],  where 
various  energies  associated  with  model  gives  rise  to  internal  and  external  forces.  The 
equilibrium  state  of  the  model  is  a  model  position  where  the  internal  deformation 
force  becomes  equal  to  the  externally  applied  force.  These  physics-based  deformable 
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models  are  useful  for  modeling  where  the  modeler  can  deform  a  surface  by  applying 

synthesized  forces,  and  for  data  fitting  where  external  forces  are  synthesized  from  a 

given  data  set  such  that  the  model  approximates  the  given  data  at  equilibrium. 

A  deformable  surface  model  is  typically  parameterized  over  the  domain  [0,  l]2. 

Let  s(u,  v,  p)  be  a  deformable  surface  model,  where  0  <  u,  v  <  1  and  p  is  a  collection 

of  the  degrees  of  freedom  p*,  i  =  1, 2, . . . ,  n  associated  with  the  model  (assuming  the 

model  has  n  degrees  of  freedom).  The  degrees  of  freedom  pi(t)  is  a  set  of  generalized 

coordinates  which  are  functions  of  time  and  are  assembled  into  the  vector  p.    Let 

fi(t)  be  the  generalized  force  represented  by  the  vector  fp  and  acting  on  p{.  Let  T  be 

the  kinetic  energy,  F  be  the  dissipation  energy  and  U  be  the  potential  energy  of  the 

deformable  surface  model.  Then,  the  Lagrangian  equation  of  motion  for  the  model 

can  be  expressed  as 

d  dT      dT      dF      8U       . 
dtopi      dpi      dpi      dpi 

Let  /x(u,  v)  be  the  mass  density  function  of  the  surface.  Then  the  kinetic  energy 
of  the  surface  is 

T    =    )-[  ( p.sTsdudv       =       |pTMp,  (2.2) 

where  M  is  called  the  mass  matrix.    Similarly,  let  7(u,  v)  be  the  damping  density 
function  of  the  surface.  Then  the  dissipation  energy  is 


F 

? 


l-  J  j  tsTsdudv      =       iprDp,  (2.3) 
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where  D  is  called  the  damping  matrix.   The  potential  energy  of  the  model  can  be 
defined  using  the  thin-plate-under-tension  energy  model  [96],  and  is  given  by 


f  ft      ?*Li*  ^—       g    d2sr  d2s 

J  J  du  du  dv  dv  d2u  d2u 

ouav  ouov  a'v  o'v  I 


where  au(u,  v)  and  0ij(u,  v)  are  elasticity  functions  which  control  tension  and  rigidity 
respectively  of  the  deformable  surface  model,  and  K  is  known  as  the  stiffness  matrix. 
The  discretized  equation  of  motion  can  be  derived  using  the  expressions  of  the  kinetic, 
dissipation  and  potential  energy  in  Eqn.2.1,  which  is  given  by 

Mp  +  Dp  +  Kp  =  f„  (2.5) 

where  fp  is  the  generalized  force  vector. 

The  free-form  deformable  models  discussed  above  were  first  introduced  to  com- 
puter graphics  and  visualization  in  Terzopoulos  et  al.  [99]  and  further  developed  by 
Terzopoulos  and  Fleischer  [97],  Pentland  and  Williams  [72],  Metaxas  and  Terzopou- 
los [65]  and  Vemuri  and  Radisavljevic  [104].  Celniker  and  Gossard  [11]  developed  a 
system  for  interactive  free-form  design  based  on  the  finite  element  optimization  of  en- 
ergy functionals  proposed  in  Terzopoulos  and  Fleischer  [97],  Bloor  and  Wilson  [8,  9], 
Celniker  and  Welch  [12]  and  Welch  and  Witkin  [105]  proposed  deformable  B-spline 
curves  and  surfaces  which  can  be  designed  by  imposing  the  shape  criteria  via  the 
minimization  of  the  energy  functionals  subject  to  hard  or  soft  geometric  constraints 
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through  Lagrange  multipliers  or  penalty  methods.  Qin  and  Terzopoulos  [77,  78,  100] 
developed  dynamic  NURBS  (D-NURBS)  which  are  very  sophisticated  models  suit- 
able for  representing  a  wide  variety  of  free-form  as  well  as  standard  analytic  shapes. 
The  D-NURBS  have  the  advantage  of  interactive  and  direct  manipulation  of  NURBS 
curves  and  surfaces,  resulting  in  physically  meaningful  hence  intuitively  predictable 
motion  and  shape  variation. 

Deformable  models  are  also  widely  used  for  shape  recovery,  segmentation,  mo- 
tion tracking  and  other  computer  vision  and  medical  imaging  applications.  A  detailed 
survey  of  deformable  models  used  in  these  techniques  can  be  found  in  Mclnerney  and 
Terzopoulos  [64]  and  the  references  therein.  Some  specific  existing  deformable  models 
used  for  shape  recovery  and  non-rigid  motion  estimation  will  be  reviewed  in  Section 
2.4. 

2.3     Shape  Modeling  Using  Physics-based  Subdivision-surface  Model 

Recursive  subdivision  surfaces  are  powerful  for  representing  smooth  geometric 
shapes  of  arbitrary  topology.  However,  they  constitute  a  purely  geometric  represen- 
tation, and  furthermore,  conventional  geometric  modeling  with  subdivision  surfaces 
may  be  difficult  for  representing  highly  complicated  objects.  For  example,  modelers 
are  faced  with  the  tedium  of  indirect  shape  modification  and  refinement  through  time- 
consuming  operations  on  a  large  number  of  (most  often  irregular)  control  vertices 
when  using  typical  subdivision  surface-based  modeling  schemes.  Despite  the  advent 
of  advanced  3D  graphics  interaction  tools,  these  indirect  geometric  operations  remain 
non-intuitive  and  laborious  in  general.  In  addition,  it  may  not  be  enough  to  obtain 
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the  most  "fair"  surface  that  interpolates  a  set  of  (ordered  or  unorganized)  data  points. 
A  certain  number  of  local  features  such  as  bulges  or  inflections  may  be  strongly  de- 
sired while  requiring  geometric  objects  to  satisfy  global  smoothness  constraints  in 
geometric  modeling  and  computer  graphics  applications.  In  contrast,  physics-based 
modeling  provides  a  superior  approach  to  shape  modeling  that  can  overcome  most  of 
the  limitations  associated  with  traditional  geometric  modeling  approaches.  Free-form 
deformable  models  governed  by  the  laws  of  continuum  mechanics  are  of  particular 
interest  in  this  context.  These  dynamic  models  respond  to  externally  applied  forces 
in  a  very  intuitive  manner.  The  dynamic  formulation  marries  the  model  geometry 
with  time,  mass,  damping  and  constraints  via  a  force  balance  equation.  Dynamic 
models  produce  smooth,  natural  motions  which  are  easy  to  control.  In  addition,  they 
facilitate  interaction  -  especially  direct  manipulation  of  complex  geometries.  Fur- 
thermore, the  equilibrium  state  of  the  model  is  characterized  by  a  minimum  of  the 
deformation  energy  of  the  model  subject  to  the  imposed  constraints.  The  deformation 
energy  functionals  can  be  formulated  to  satisfy  local  and  global  modeling  criteria,  and 
geometric  constraints  relevant  to  shape  design  can  also  be  imposed.  The  dynamic 
approach  subsumes  all  of  the  aforementioned  modeling  capabilities  in  a  formulation 
which  grounds  everything  in  real-world  physical  behavior. 

A  severe  limitation  of  the  existing  deformable  models,  including  D-NURBS,  is 
that  they  are  defined  on  a  rectangular  parametric  domain.  Hence,  it  can  be  very 
difficult  to  model  surfaces  of  arbitrary  genus  using  these  models.  In  a  recent  work, 
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DeRose  et  al.  [25]  assigned  material  properties  to  control  meshes  at  various  subdivi- 
sion levels  in  order  to  simulate  cloth  dynamics  using  subdivision  surfaces.  Note  that, 
they  assign  physical  properties  on  the  control  meshes  at  various  levels  of  subdivision 
and  not  on  the  limit  surface  itself,  and  hence  can  not  solve  the  modeling  goal  we 
are  trying  to  achieve.  In  this  dissertation,  a  dynamic  framework  is  presented  for 
subdivision  surfaces  which  combines  the  benefits  of  subdivision  surfaces  for  modeling 
arbitrary  topology  as  well  as  that  of  dynamic  splines  for  interactive  shape  manip- 
ulation by  applying  synthesized  forces.  The  proposed  dynamic  framework  presents 
the  modeler  a  formal  mechanism  of  direct  and  intuitive  manipulation  of  the  smooth 
limit  surface,  as  if  they  were  seamlessly  sculpting  a  piece  of  real-world  "clay."  A 
dynamic  framework  for  subdivision  surface-based  wavelets  adds  flexibility  in  the  pro- 
posed modeling  paradigm.  It  enables  a  multiresolution  representation  of  the  evolving 
control  mesh  defining  the  smooth  limit  surface,  and  the  modeler  can  control  the 
extent  of  manipulation  effect  by  choosing  a  proper  editing  level.  The  formulation 
and  implementation  details  of  these  dynamic  frameworks  are  discussed  in  subsequent 
chapters. 
2.4    Shape  and  Motion  Estimation  Using  Phvsics-based  Subdivision-surface  Model 

The  dynamic  subdivision  surface  model  has  been  developed  primarily  for  mod- 
eling arbitrary  (known)  topology  where  modelers  can  directly  manipulate  the  limit 
surface  by  applying  synthesized  forces  in  an  intuitive  fashion.  However,  another  im- 
portant application  of  the  dynamic  subdivision  surfaces  is  in  non-rigid  shape  and 
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motion  reconstruction/recovery.  Accurate  shape  recovery  requires  distributed  pa- 
rameter models  which  typically  possess  a  large  number  of  degrees  of  freedom.  On 
the  other  hand,  efficient  shape  representation  imposes  the  requirement  of  geometry 
compression,  i.e.,  models  with  fewer  degrees  of  freedom.  These  requirements  are 
conflicting  and  numerous  researchers  have  been  seeking  to  strike  a  balance  between 
these  contradicting  requirements  [5,  21,  47,  49,  62,  89,  100,  104].  Another  important 
criterion  in  the  design  of  shape  models  is  that  the  initialization  of  the  model  during 
the  shape  recovery  process  should  not  be  restricted  to  globally  parameterized  input 
meshes  since  it  may  be  infeasible  to  globally  parameterize  shapes  of  arbitrary  topol- 
ogy. A  physics-based  model  best  satisfying  the  above  mentioned  criteria  is  an  ideal 
candidate  for  a  solution  to  the  shape  recovery  problem  for  obvious  reasons. 

Deformable  models  which  come  in  many  varieties,  have  been  used  to  solve  this 
problem  in  the  physics-based  modeling  paradigm.  These  models  involve  the  use  of 
either  fixed  size  [21,  66,  71,  98,  104]  or  adaptive  size  [15,  41,  47,  62,  95,  101]  grids. 
The  models  with  fixed  grid  size  generally  use  less  number  of  degrees  of  freedom  for 
representation  at  the  cost  of  accuracy  of  the  recovered  shape.  On  the  other  hand, 
models  using  adaptive  grids  generally  need  large  number  of  degrees  of  freedom  to 
recover  the  shapes  accurately.  Note  that  the  shapes  being  recovered  from  the  image 
data  are  smooth  in  most  of  the  medical  applications,  i.e.  the  anatomical  structures 
even  with  considerable  amount  of  details  have  more  or  less  a  C1  surface.  Under 
these  circumstances,  the  finite  element  approaches  as  in  Cohen  and  Cohen  [21]  and 
Mclnerney  and  Terzopoulos  [62]  need  a  /arge  number  of  degrees  of  freedom  for  deriving 
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a  smooth  and  accurate  representation.  In  addition,  they  can  not  represent  shapes 
with  known  arbitrary  topology.  Moreover,  almost  all  of  these  schemes  require  a 
globally  parameterized  mesh  as  their  input  which  may  be  infeasible  at  times. 

The  proposed  model  solves  the  shape  recovery  problem  very  effectively  as  it 
can  recover  shapes  from  large  range  and  volume  data  sets  using  very  few  degrees 
of  freedom  (control  vertices)  for  its  representation  and  can  cope  with  any  arbitrary 
input  mesh,  not  necessarily  parameterized.  The  initialized  model  deforms  under  the 
influence  of  synthesized  forces  to  fit  the  data  set  by  minimizing  its  energy.  Once  the 
approximate  shape  is  recovered,  the  model  is  further  subdivided  automatically  and  a 
better  approximation  to  the  input  data  set  is  achieved  using  more  degrees  of  freedom. 
The  process  of  subdivision  after  achieving  an  approximate  fit  is  continued  till  a  pre- 
scribed error  criteria  for  fitting  the  data  points  is  achieved.  The  proposed  dynamic 
subdivision  surface  models  have  also  been  successfully  used  in  motion  tracking  and 
visualization  of  a  dynamically  deforming  shape  from  a  time  sequence  of  volumetric 
data  sets. 

2.5     Unified  Framework  for  Shape  Recovery  and  Shape  Modeling 

Currently,  shape  recovery  and  shape  modeling  are  viewed  as  two  distinct  areas  in 
computer  graphics  and  geometric  modeling  literature.  However,  there  are  potential 
benefits  if  these  two  can  be  combined  in  an  unified  framework.  For  example,  the 
modeler  starts  from  scratch  to  build  a  specific  model  in  a  typical  geometric  modeling 
scenario.  First,  a  rough  shape  is  modeled  and  then  it  is  fine-tuned  by  manipulating 
control  vertex  positions  to  obtain  the  desired  effects.  This  turns  out  a  cumbersome 


25 


process  in  general.  On  the  other  hand,  shape  recovery  using  state-of-the-art  methods 
yield  large  polygonal  meshes  which  are  very  difficult  to  manipulate,  especially  for 
global  changes  in  shape.  Most  often  the  resulting  meshes  from  a  shape  recovery 
application  are  not  directly  amenable  to  multiresolution  analysis.  Computationally 
expensive  re-meshing  techniques  are  needed  to  convert  these  meshes  into  a  specific 
type  of  meshes  on  which  multiresolution  analysis  can  be  performed.  This  specific  type 
of  mesh  is  known  as  mesh  with  "subdivision-connectivity",  implying  a  topologically 
equivalent  mesh  with  the  same  connectivity  as  of  the  given  mesh  can  be  obtained  by 
recursive  subdivision  of  a  very  simple  known  initial  mesh. 

The  proposed  dynamic  framework  combines  shape  recovery  and  shape  modeling 
in  an  unified  framework  where  the  modeler  can  scan  in  3D  points  of  a  prototype 
model,  recover  the  shape  using  the  proposed  dynamic  subdivision  surface  model,  and 
edit  at  any  desired  resolution  using  physics-based  force  tools.  Thus,  the  modeler  is  re- 
lieved of  the  burden  of  both  building  an  initial  model  and  editing  a  cumbersome  huge 
polygonal  mesh.  The  shape  recovery  process  starts  with  a  smooth  subdivision  surface 
model  which  has  a  simple  initial  mesh.  This  model  is  deformed  using  synthesized 
forces  from  the  given  data  points  and  is  automatically  refined  using  some  pre-defined 
error  in  fit  criteria.  The  initial  mesh  of  the  final  recovered  smooth  shape  has  subdivi- 
sion connectivity  in-built,  as  it  is  obtained  by  subdivision  refinement,  and  therefore 
no  re-meshing  is  needed  for  multiresolution  analysis.  The  dynamic  framework  for 
subdivision  surface-based  wavelets  makes  multiresolution  editing  using  physics-based 
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force  tools  easy  to  perform.  These  advantages  of  shape  modeling  and  shape  recovery 
in  a  unified  framework  will  be  further  illustrated  in  later  chapters. 


CHAPTER  3 
DYNAMIC  CATMULL-CLARK  SUBDIVISION  SURFACES 


Subdivision  surfaces,  as  mentioned  earlier,  can  be  broadly  classified  into  two 
categories  -  approximating  and  interpolatory  subdivision  schemes.  The  approxi- 
mating subdivision  schemes  typically  generalize  recursive  subdivision  techniques  for 
generating  limit  surfaces  with  known  parameterizations.  In  this  chapter,  a  dynamic 
framework  will  be  presented  for  Catmull-Clark  subdivision  scheme,  a  popular  approx- 
imating subdivision  technique.  The  Catmull-Clark  subdivision  scheme  generalizes  the 
idea  of  obtaining  a  bicubic  B-spline  surface  by  recursive  subdivision  of  a  rectangu- 
lar initial  mesh.  Before  discussing  the  dynamic  framework,  first  the  Catmull-Clark 
subdivision  scheme  is  briefly  reviewed. 

3.1     Overview  of  the  Catmull-Clark  Subdivision  Scheme 

Catmull-Clark  subdivision  scheme,  like  any  other  subdivision  scheme,  starts 
with  a  user-defined  mesh  of  arbitrary  topology.  It  refines  the  initial  mesh  by  adding 
new  vertices,  edges  and  faces  with  each  step  of  subdivision  following  a  fixed  set  of 
subdivision  rules.  In  the  limit,  a  recursively  refined  polyhedral  mesh  will  converge  to 
a  smooth  surface.  The  subdivision  rules  are  as  follows: 
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(a) 


(b) 


(c) 


Figure  3.1.  Catmull-Clark  subdivision  :  (a)  initial  mesh,  (b)  mesh  obtained  after  one 
step  of  Catmull-Clark  subdivision,  and  (c)  mesh  obtained  after  another  subdivision 
step. 


•  For  each  face,  a  new  face  vertex  is  introduced  which  is  the  average  of  all  the 
old  vertices  defining  the  face. 

•  For  each  (non-boundary)  edge,  a  new  edge  vertex  is  introduced  which  is  the 
average  of  the  following  four  points:  two  old  vertices  defining  the  edge  and  two 
new  face  vertices  of  the  faces  adjacent  to  the  edge. 

•  For  each  (non-boundary)  vertex  V,  a  new  vertex  is  introduced  whose  position 
is  £  +  2E  +  \n~  >  ,  where  F  is  the  average  of  the  new  face  vertices  of  all  faces 
adjacent  to  the  old  vertex  V,  E  is  the  average  of  the  midpoints  of  all  the  edges 
incident  on  the  old  vertex  V  and  n  is  the  number  of  the  edges  incident  on  the 
vertex. 

•  New  edges  are  formed  by  connecting  each  new  face  vertex  to  the  new  edge 
vertices  of  the  edges  defining  the  old  face  and  by  connecting  each  new  vertex 
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corresponding  to  an  old  vertex  to  the  new  edge  vertices  of  all  old  edges  incident 
on  the  old  vertex. 

•  New  faces  are  defined  as  faces  enclosed  by  new  edges. 

An  example  of  Catmull-Clark  subdivision  on  an  initial  mesh  is  shown  in  Fig.3.1. 
The  most  important  property  of  the  Catmull-Clark  subdivision  surfaces  is  that  a 
smooth  surface  can  be  generated  from  any  control  mesh  of  arbitrary  topology.  There- 
fore, this  subdivision  scheme  is  extremely  valuable  for  modeling  various  complicated 
geometric  objects  of  arbitrary  topology.  Catmull-Clark  subdivision  surfaces  include 
standard  bicubic  B-spline  surfaces  as  their  special  case  (i.e.,  the  limit  surface  is  a 
bicubic  B-spline  surface  for  a  rectangular  mesh  with  all  non-boundary  vertices  of 
degree  4).  In  addition,  the  aforementioned  subdivision  rules  generalize  the  recur- 
sive bicubic  B-spline  patch  subdivision  algorithm.  For  non-rectangular  meshes,  the 
limit  surface  converges  to  a  bicubic  B-spline  surface  except  at  a  finite  number  of  ex- 
traordinary points.  These  extraordinary  points  correspond  to  extraordinary  vertices 
(vertices  whose  degree  is  not  equal  to  4)  in  the  mesh.  Note  that,  after  the  first  sub- 
division, all  faces  are  quadrilaterals,  hence  all  the  new  vertices  created  subsequently 
will  have  four  incident  edges.  The  number  of  extraordinary  points  on  the  limit  sur- 
face is  a  constant,  and  is  equal  to  the  number  of  extraordinary  vertices  in  the  refined 
mesh  obtained  after  applying  one  step  of  the  Catmull-Clark  subdivision  on  the  initial 
mesh.  The  limit  surface  is  curvature-continuous  everywhere  except  at  extraordinary 
vertices,  where  only  tangent  plane  continuity  is  achieved.  In  spite  of  the  popularity 
of  Catmull-Clark  subdivision  surfaces  for  representing  complex  geometric  shapes  of 
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arbitrary  topology,  these  subdivision  surfaces  may  not  be  easily  parameterizable  and 
deriving  a  closed-form  analytic  formulation  of  the  limit  surface  can  be  very  difficult. 
These  deficiencies  preclude  their  immediate  pointwise  manipulation  and  hence  may 
restrain  the  applicability  of  these  schemes.  In  the  rest  of  the  chapter,  a  dynamic 
framework  is  developed  for  the  Catmull-Clark  subdivision  surfaces  which  offers  a 
closed-form  analytic  formulation  and  allows  users  to  manipulate  the  model  directly 
and  intuitively.  It  may  be  noted  that  recently  a  scheme  was  proposed  by  Stam  [90] 
for  directly  evaluating  the  Catmull-Clark  subdivision  surfaces  at  arbitrary  parameter 
values,  but  the  work  presented  here  on  dynamic  Catmull-Clark  subdivision  surfaces 
[57,  60,  76]  was  published  prior  to  the  publication  of  the  above-mentioned  work. 

3.2     Formulation 

In  this  section,  a  systematic  formulation  of  the  new  dynamic  model  based  on 
Catmull-Clark  subdivision  surfaces  is  presented.  This  dynamic  model  treats  the 
smooth  limit  surface  as  a  function  of  its  initial  mesh.  However,  the  control  vertex  po- 
sitions need  to  be  updated  continually  at  every  level  of  subdivision  in  order  to  develop 
the  dynamic  framework.  Note  that,  all  the  vertices  introduced  through  subdivision 
are  obtained  as  an  affine  combination  of  control  vertex  positions  in  the  initial  mesh. 
Therefore,  the  dynamic  behavior  of  the  limit  surface  can  be  controlled  by  formulat- 
ing the  dynamic  model  on  the  initial  mesh  itself,  the  only  exception  being  the  case 
when  the  initial  mesh  has  non-rectangular  faces.  This  problem  can  be  circumvented 
by  taking  the  mesh  obtained  through  one  step  of  subdivision  as  the  initial  mesh. 
It  may  be  noted  that  an  additional  subdivision  step  may  be  needed  in  some  cases 
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to  isolate  the  extraordinary  points,  and  the  resulting  mesh  is  treated  as  the  initial 
mesh.  A  typical  example  of  the  above  mentioned  scenario  is  when  the  initial  mesh  is 
a  tetrahedron  where  two  subdivision  steps  are  needed,  and  the  dynamic  framework 
can  be  developed  by  treating  the  mesh  obtained  after  two  subdivision  steps  as  the 
initial  mesh.  To  define  the  limit  surface  using  the  vertices  of  the  initial  mesh,  the 
enumeration  of  the  bicubic  patches  in  the  limit  surface  is  necessary.  In  the  next  two 
sections,  various  schemes  are  presented  for  assigning  the  bicubic  patches  to  various 
faces  of  the  initial  mesh. 

3.2.1 Assigning  Patches  to  Regular  Faces 

In  Fig.3.2,  a  rectangular  control  mesh  is  shown  along  with  the  bicubic  B-spline 
surface  (4  patches)  in  the  limit  after  an  infinite  number  of  subdivision  steps.  Note 
that,  each  of  the  bicubic  patches  in  the  limit  surface  is  defined  by  a  rectangular  face 
with  each  vertex  of  degree  four,  thereby  accounting  for  16  control  points  (from  its 
8  connected  neighborhood)  needed  to  define  a  bicubic  surface  patch  in  the  limit. 
Therefore,  for  each  rectangular  face  in  the  initial  mesh  with  a  degree  of  4  at  each 
vertex,  the  corresponding  bicubic  surface  patch  can  be  assigned  to  it  in  a  straight 
forward  way.  In  Fig.3.2,  the  surface  patches  Si,S2,S3  and  54  are  assigned  to  face 
Fi,  F2,  F3  and  F4  respectively.  The  16  control  points  for  the  patch  Si,  corresponding 
to  face  F, ,  are  highlighted  in  Fig.3.2.  Note  that,  the  initial  control  mesh  can  be  viewed 
as  the  parametric  domain  of  the  limit  surface.  Therefore,  face  Fi  can  be  thought  of 
as  the  portion  of  the  parametric  domain  over  which  the  patch  Si  is  defined,  i.e.,  has 
non-zero  values.  Nevertheless,  each  rectangular  face  (e.g.  F%)  can  be  parametrically 
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Figure  3.2.  A  rectangular  mesh  and  its  limit  surface  consisting  of  4  bicubic  surface 
patches. 
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defined  over  [0, 1]  ,  and  hence,  all  bicubic  B-spline  patches  defined  by  16  control 
points  are  locally  parameterized  over  [0, 1]  . 

3.2.2     Assigning  Patches  to  Irregular  Faces 

In  Fig.3.3,  a  mesh  containing  an  extraordinary  point  of  degree  3  and  its  limit  sur- 
face are  shown.  The  faces  Fo,  Fi,  . . . , Fg  are  assigned  to  bicubic  patches  So,  Si,  ■  ■  ■ ,  Sg 
respectively  (as  they  all  have  vertices  of  degree  4)  following  the  aforementioned 

scheme.  However,  the  central  smooth  surface  enclosed  by  the  patches  So,  Sj S$ 

consists  of  infinite  number  of  bicubic  patches  converging  to  a  point  in  the  limit.  A  re- 
cursive way  of  enumerating  these  bicubic  patches  and  assigning  them  to  various  faces 
at  different  levels  need  to  be  developed  in  order  to  formulate  the  dynamic  framework 
for  Catmull-Clark  subdivision  surface  model. 

The  idea  of  enumerating  the  bicubic  patches  corresponding  to  faces  having  an 
extraordinary  vertex  is  shown  in  Fig.3.4  where  a  local  subdivision  of  the  mesh  con- 
sisting of  faces  Fo,  F\, ...,  F8,  P0,  Pi,P?  (and  not  the  other  boundary  faces)  of  Fig.3.3 
is  carried  out.  Topological^,  the  resulting  local  subdivision  mesh  (shown  as  dotted 
mesh)  is  exactly  the  same  as  the  mesh  in  Fig.3.3  and  hence  exactly  the  same  number 
of  bicubic  patches  can  be  assigned  to  its  faces  with  vertices  of  degree  4  as  is  evident 
from  Fig.3.4  (the  new  faces  and  the  corresponding  patches  are  marked  by  "p"  and 
"n"  respectively).  This  process  of  local  subdivision  and  assignment  of  bicubic  patches 
around  an  extraordinary  point  can  be  carried  out  recursively  and  in  the  limit,  the 
enclosed  patch  corresponding  to  faces  sharing  the  extraordinary  point  will  converge 
to  a  point.  However,  there  is  no  need  to  carry  out  an  infinite  number  of  subdivision 
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Figure  3.3.  A  mesh  with  an  extraordinary  point  of  degree  3  and  its  limit  surface. 
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Figure  3.4.  Local  subdivision  around  the  extraordinary  point  and  the  limit  surface. 
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steps.  This  description  is  for  formulation  purposes  only  and  the  exact  implementation 
will  be  detailed  in  a  later  section. 

3.2.3     Kinematics 

In  this  section,  the  mathematics  for  the  kinematics  of  the  limit  surface  is  devel- 
oped via  illustrative  examples  and  then  the  generalized  formulas  are  presented.  It 
may  be  noted  that  a  single  bicubic  B-spline  patch  is  obtained  as  the  limiting  process 
of  the  Catmull-Clark  subdivision  algorithm  applied  to  an  initial  4  by  4  rectangular 
control  mesh.  Let  Sp(u,v),  where  (u,  v)  e  [0,  l]2,  denote  this  bicubic  B-spline  patch 
which  can  be  expressed  analytically  as 


sp(u,u)  =  (x(u,v),y(u,v),z(u,v)f  =  J2YlA'JBiAu)BiAv)'  t3-1) 

i=0  j=0 


where  dM  represents  a  3-dimensional  position  vector  at  the  (i,j)th  control  point 
location  and  Bhi(u)  and  Bjti(v)  are  the  cubic  B-spline  basis  functions.  The  subscript 
p  on  s  denotes  the  patch  under  consideration.  Expressing  Eqn.3.1  in  a  generalized 
coordinate  system  we  have 

sp  =  Jpqp,  (3.2) 

where  3P  is  a  transformation  matrix  storing  the  basis  functions  of  a  bicubic  B-spline 
patch,  and  is  of  size  (3,48).  Vector  qp  is  the  concatenation  of  all  control  points 
defining  a  B-spline  patch  in  3D.  Note  that  in  the  concatenation  of  the  control  points, 
each  control  point  has  an  (x,y,  z)  component.  For  example,  the  (x,y,z)  components 
of  the  control  point  (i,j)  correspond  to  positions  3£,3fc  +  l,3fc  +  2  -  where,  k  = 
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4i  +  j  -  respectively  in  the  vector  (jp.  We  can  express  the  entries  of  3P  explicitly 
in  the  following  way:  Jp(0,  k)  =  Jp(l,fc  +  1)  =  Jp(2, fc  +  2)  =  Biti{u)Bjti{v)  and 
Jp(0,  k  +  1)  =  Jp(0,  k  +  2)  =  Jp(l,  fc)  =  Jp(l,fc  +  2)  =  Jp(2,  k)  =  Jp(2,  k  +  1)  =  0. 

Limit  surface  with  many  bicubic  patches  from  a  rectangular  initial  mesh 

Now  a  limit  surface  consisting  of  many  bicubic  surface  patches  obtained  after 
applying  an  infinite  number  of  subdivision  steps  to  a  rectangular  initial  mesh  is 
considered.  For  example,  let  the  limit  surface  of  Fig.3.2  be  sm  which  can  be  written 


sm(H,u)  =  STni(2u,2w)  +  sm2(2(u--),2!;)  +  sm3(2(«--),2();--))+sm4(2«,2(i;--)), 

(3.3) 
where  smi(2u,2v)  =  sm(u,v)  for  0  <  u,v  <  5,  and  0  otherwise.  Similarly,  sm2,sms 
and  smi  are  also  equal  to  sm(u,v)  for  an  appropriate  range  of  values  of  u,v  and  0 
outside.  It  may  be  noted  that  smi,sm2,sm3,sra4  correspond  to  patches  Sl,S2,S3,Si 
respectively  in  Fig.3.2.  Eqn.3.3  in  generalized  coordinates  becomes 


sra  =  Jiqi  +  J2q2  +  J3Q3  +  J4q4  =  XI  J'*'  (3-4) 


4 

<=1 


where  J^s  are  the  transformation  matrices  of  size  (3,48)  and  qjS  are  the  (x,y,z)  com- 
ponent concatenation  of  a  subset  of  the  control  points  of  sm  defining  smi,  i  =  1,2,3 


38 


and  4.  A  more  general  expression  for  sm  is 

4 

sm  =  Ji  A[qm  +  J2A2qm  +  J3  A3qm  +  J4  A.,qm  =  £  J,  A;qm  =  Jraqm.         (3.5) 

i=i 

Where,  qm  is  the  75-component  vector  of  3D  positions  of  the  25  vertex  control  mesh 
defining  the  limit  surface  sm.    Matrices  At,  1  <  i  <  4,  are  of  size  (48, 75)  ,  each 
row  consisting  of  a  single  nonzero  entry  (=  1)  and  the  (3,  75)-sized  matrix  JTO  = 
Em  J.A,. 
Limit  surface  with  many  bicubic  patches  from  an  arbitrary  initial  mesh 

The  stage  is  now  set  to  define  the  limit  surface  s  using  the  vertices  of  initial 
mesh  M  for  any  arbitrary  topology,  assuming  all  faces  are  quadrilateral  and  no  face 
contains  more  than  one  extraordinary  point  as  its  vertex  (i.e.,  extraordinary  points 
are  isolated).  As  mentioned  earlier,  if  these  assumptions  are  not  satisfied,  one  or  two 
steps  of  global  subdivision  may  be  required  and  the  resulting  mesh  can  be  treated 
as  the  initial  mesh.  Let  the  number  of  vertices  in  the  initial  mesh  M  be  a,  and  let 
/  of  these  be  the  extraordinary  vertices.  The  number  of  faces  in  the  initial  mesh  are 
assumed  to  be  b,  and  k  of  these  faces  have  vertices  with  degree  4  (henceforth  termed  a 
"normal  face")  and  each  of  the  remaining  (b—k)  faces  have  one  of  the  /  extraordinary 
vertices  (henceforth  termed  a  "special  face").  Let  p  be  the  3a  =  N  dimensional  vector 
containing  the  control  vertex  positions  in  3D.  Using  the  formulations  in  Section  3.2.1 
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and  Section  3.2.2,  the  smooth  limit  surface  can  be  expressed  as 


■=2>+E«».  (3-6) 

i=l  j=l 


where  n^  is  a  single  bicubic  patch  assigned  to  each  of  the  normal  faces  and  Sj  is  a 
collection  of  infinite  number  of  bicubic  patches  corresponding  to  each  of  the  extraor- 
dinary points.  As  s  is  a  surface  defined  over  the  faces  of  the  initial  control  mesh,  each 
11^  can  be  viewed  as  a  bicubic  patch  defined  over  the  corresponding  regular  face  in  the 
initial  control  mesh.  Similarly,  each  Sj  can  be  defined  over  the  corresponding  irregu- 
lar faces  in  the  initial  mesh  (refer  to  Fig.3.2  -  Fig.3.4  for  the  detailed  description  on 
parametric  domains  of  subdivision  surfaces).  For  the  simplicity  of  the  following  math- 
ematical derivations  on  the  dynamic  model,  from  now  on  the  parametric  domain  will 
not  be  explicitly  provided  in  the  formulation.  Employing  the  same  approach  taken 
before  to  derive  Eqn.3.5,  it  can  be  shown  that 


£  n<  =  ETOCp,)  =  (£C  J.)("A,))P  =  ("J)p,  (3.7) 

i=i  t=i  i=i 


where  nJi,np,  and  "A,  are  the  equivalent  of  J,,p,  in  Eqn.3.4  and  A,  in  Eqn.3.5  respec- 
tively. The  pre-superscript  n  is  used  to  indicate  that  these  mathematical  quantities 
describe  bicubic  patch  in  the  limit  surface  corresponding  to  normal  faces. 

The  following  notational  convention  will  be  used  for  describing  various  math- 
ematical quantities  used  in  the  derivation  of  the  expression  for  a  collection  of  infinite 
number  of  bicubic  patches  around  an  extraordinary  vertex.    The  pre-superscript  s 


40 

is  used  to  represent  a  collection  of  bicubic  patches  around  an  extraordinary  vertex, 
the  subscript  j  is  used  to  indicate  the  j-th  extraordinary  point,  the  post-superscript 
represents  the  exponent  of  a  mathematical  quantity  and  the  level  indicator  (to  rep- 
resent various  levels  of  subdivision  in  the  local  control  mesh  around  an  extraordinary 
vertex)  is  depicted  via  a  subscript  on  the  curly  braces. 

The  expression  for  Sj  is  derived  using  the  recursive  nature  of  local  subdivision 
around  an  extraordinary  vertex  as  shown  in  Section  3.2.2.  First,  Sj  can  be  expressed 
as 

where  the  first  term  of  Eqn.3.8  is  the  generalized  coordinate  representation  of  the 
bicubic  B-spline  patches  corresponding  to  the  normal  faces  of  the  new  local  subdivi- 
sion mesh  obtained  after  one  subdivision  step  on  the  local  control  mesh  (similar  to 
those  patches  marked  n  in  Fig. 3. 4).  {Sj}l  represents  the  rest  of  the  infinite  bicubic 
B-spline  patches  surrounding  the  extraordinary  point  (similar  to  the  central  patch 
enclosed  by  patches  marked  n  in  Fig.3.4).  The  vertices  in  the  newly  obtained  lo- 
cal subdivision  mesh  {*Pj}j  can  be  expressed  as  a  linear  combination  of  a  subset 
of  the  vertices  of  the  initial  mesh  M.  (which  will  contribute  to  the  local  subdivi- 
sion) following  the  subdivision  rules.  We  can  name  this  subset  of  initial  control 
vertices  {"Pj}0-  Furthermore,  there  exists  a  matrix  {sBj}j  of  size  (3c,  3d),  such  that 
{^BjJjCpjJij  =  {sPj},  where  {spJ}1  and  {sPj}0  are  vectors  of  dimension  3c  and  3d 
respectively.  Applying  the  idea  of  recursive  local  subdivision  again  on  {sj},,  Sj  can 
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be  further  expanded  as 

Sj  =  {'J,}l{'Bi}1{*Pi}9'+{,J,}2{*%}a{'p,}]  +  (%}r  (3.9) 

In  the  above  derivation,  {"Pj}!  is  a  vector  of  dimension  3d,  comprising  of  a 
subset  of  the  vertices  defining  the  3c  dimensional  vector  {"Pj}1-  Note  that,  {"Pj}, 
has  the  same  structure  as  {sp:)}0,  therefore,  there  exists  a  (3d, 3d)  matrix  {"Cj}l 
such  that  {*Ci}j{*Pj},|  =  {aPj}r  Each  subdivision  of  a  local  mesh  with  d  vertices 
creates  a  new  local  mesh  with  c  vertices  which  contributes  a  fixed  number  of  bicubic 
B-spline  patches.  Proceeding  one  step  further,  it  can  be  shown  that 

h  =  f'JArBACpA  +  C^U'BiU'c^fp,), 

+{°J]}3{'Bj}3{sCj}21{sp}}0  +  {■/},.  (3.10) 

Because  of  the  intrinsic  property  of  the  local  recursive  subdivision  around  the  ex- 
traordinary point,  we  have  {"3j}l  =  {sJj}2  =  •■■  =  {sJj}„  =  ...  =  {'3j}x-  In 
addition,  the  subdivision  rules  remain  the  same  throughout  the  refinement  process, 
we  also  have  {'B^j  =  {sBj}2  =  ...  =  {sBj}n  =  ...  =  {'B^.  So,  the  above 
equations  can  be  further  simplified  leading  to 

s,    =    {,JJ}1{sBJ}1{'pJ}0+{"JJ}1{"BJ}l{»CJ}1{spj}0 
+{sJJ}1{sBJ}1{sCJ}*{>pJ}0  +  ... 

oo 

=    {*Ji}l{*Bi}1(£{*Ci}'1){'p,}0.  (3.11) 

i=0 
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Vector  Sj  can  be  written  as 

»j  =  ChWPi),  (3-12) 

where  "Jj  =  {8Jj}1{sB:,}1(^Ic^0  {*C/},)  and  *Pj  =  {sp3}0-  The  idea  of  local  recursive 
subdivision  around  an  extraordinary  point  is  illustrated  in  Fig.3.5.  Note  that,  each 
vertex  position  in  the  subdivided  mesh  is  obtained  by  an  affine  combination  of  some 
vertices  in  the  previous  level  and  hence  any  row  of  {sCj}j  sums  to  1.  The  largest 
eigenvalue  of  such  a  matrix  is  1  and  it  can  be  shown  that  the  corresponding  infinite 
series  is  convergent  following  a  similar  approach  as  in  Halstead  et  al.  [38].  The  rest 
of  the  derivation  leading  to  an  expression  for  s  is  relatively  straight  forward.  Using 
the  same  approach  used  to  derive  the  Eqn.3.7,  it  can  be  shown  that 


i:  h = Ecm'Pi)  =  (E(*j*)('Ai))P = (sj)P.        (3.13) 


From  Eqn.3.6,  Eqn.3.7  and  Eqn.3.13, 


s  =  ("J)p  +  ('iJ)p.  (3.14) 


Let  J  =  (nJ)  +  («J),  and  hence 


s  =  Jp.  (3.15) 
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Figure  3.5.  Local  subdivision  around  the  extraordinary  point  and  the  corresponding 
patches  in  the  limit  surface  from  different  levels  of  subdivision. 
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3.2.4     Dynamics 

In  the  previous  section,  an  expression  is  derived  for  the  smooth  limit  surface 
obtained  via  infinite  subdivision  steps.  This  expression  treats  the  smooth  limit  surface 
as  a  function  of  the  control  vertex  positions  in  the  initial  mesh.  Now  the  vertex 
positions  in  the  initial  mesh  defining  the  smooth  limit  surface  s  are  treated  as  a 
function  of  time  in  order  to  embed  the  Catmull-Clark  subdivision  model  in  a  dynamic 
framework.  The  velocity  of  this  surface  model  can  be  expressed  as 

s(u,w,p)  =  Jp,  (3.16) 

where  an  overstruck  dot  denotes  a  time  derivative.  The  physics  of  the  dynamic 
subdivision  surface  model  is  based  on  the  work-energy  version  of  Lagrangian  dynamics 
[36]  and  is  formulated  in  an  analogous  way  to  the  dynamic  framework  presented  in 
Section  2.2. 

Let  fi{u,  v)  be  the  mass  density  function  of  the  surface.  Then  the  kinetic  energy 
of  the  surface  is 

T    =    -f  J  iisrsdudv      =      -pTMp,  (3.17) 

where  (using  Eqn.3.16) 

M=  /  /  n3TUudv  (3.18) 
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is  the  mass  matrix  of  size  (N,N).    Similarly,  let  j(u,v)  be  the  damping  density 
function  of  the  surface.  Then  the  dissipation  energy  is 

F    =    -  f  J  jsTsdudv      =        2PTDp,  (3.19) 


where 


D=  /  f^JTJdudv  (3.20) 


is  the  damping  matrix  of  size  (N,  N). 

The  potential  energy  of  the  smooth  limit  surface  can  be  expressed  as 


V    =    ipTKp,  (3.21) 


where  K  is  the  stiffness  matrix  of  size  (N,N).  A  thin-plate-under-tension  energy 
model  [96]  is  used  to  compute  the  elastic  potential  energy  of  the  dynamic  subdivision 
surface.  The  corresponding  expression  for  the  stiffness  matrix  K  is 

K  =  J  j(auiTJu  +  a223l3v  +  PnJljuu  +  A*J£,J«.  +  0223lv3vv)dudv.     (3.22) 

where  the  subscripts  on  J  denote  the  parametric  partial  derivatives.  The  aa(u,v) 
and  /3ij(u,v)s  are  elasticity  functions  controlling  local  tension  and  rigidity  in  the  two 
parametric  coordinate  directions.  Note  that  the  thin-plate  energy  expression  might 
diverge  at  extraordinary  points  on  the  limit  surface  for  Catmull-Clark  subdivision 
scheme  as  shown  in  Halstead  et  al.    [38].   Several  methods  have  been  suggested  in 
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Halstead  et  al.  [38]  to  overcome  the  problem  of  the  divergent  series.  However,  in  this 
dissertation,  the  problem  is  circumvented  by  setting  the  corresponding  rigidity  coef- 
ficients to  be  zero  at  these  points.  Therefore,  the  thin-plate  energy  at  extraordinary 
points  is  zero.  The  effect  is  negligible  to  the  overall  thin-plate  energy  as  very  few 
extraordinary  points  are  present  in  the  smooth  limit  surface. 

Using  the  expressions  for  the  kinetic,  dissipation  and  potential  energy  in  Eqn.2.1, 
the  same  motion  equation  as  in  Eqn.2.5  can  be  obtained.  The  generalized  force  vector 
fp,  which  can  be  obtained  through  the  principle  of  virtual  work  [36],  is  expressed  as 

fP  =  J  J  3Tf(u,  v,  t)dudv  (3.23) 

Different  types  of  forces  can  be  applied  on  the  smooth  limit  surface,  and  the  limit 
surface  would  evolve  over  time  according  to  Eqn.2.5  to  obtain  an  equilibrium  position 
characterized  by  a  minimum  of  the  total  model  energy. 

3.2.5     Multilevel  Dynamics 

At  this  point,  a  dynamic  framework  is  developed  for  the  Catmull-Clark  subdi- 
vision scheme  where  the  smooth  limit  surface  evolves  over  time  in  response  to  the 
applied  forces.  The  entire  process  can  be  described  as  follows.  Given  an  initial  mesh, 
a  smooth  surface  is  obtained  in  the  limit.  Users  can  directly  apply  synthesized  forces 
to  this  smooth  limit  surface  to  enforce  various  functional  and  aesthetic  constraints. 
This  direct  manipulation  is  then  transferred  back  as  virtual  forces  acting  on  the  initial 
mesh  through  a  transformation  matrix  (Eqn.3.23),  and  the  initial  mesh  (as  well  as 
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the  underlying  smooth  limit  surface)  deforms  continuously  over  time  until  an  equilib- 
rium position  is  obtained.  The  deformation  of  the  surface  in  response  to  the  applied 
forces  is  governed  by  the  motion  equation  (Eqn.2.5).  Within  the  physics-based  mod- 
eling framework,  the  limit  surface  evolves  as  a  consequence  of  the  evolution  of  the 
initial  mesh.  One  can  apply  various  types  of  forces  on  the  limit  surface  to  obtain 
a  desired  effect,  but  the  possible  level  of  details  appearing  in  a  shape  that  can  be 
obtained  through  evolution  is  constrained  by  the  number  of  control  vertices  in  the 
initial  mesh.  It  might  be  necessary  to  increase  the  number  of  control  vertices  in  the 
initial  mesh  in  order  to  obtain  a  detailed  shape  through  this  evolution  process. 

The  number  of  control  vertices  defining  the  same  smooth  limit  surface  can  be 
increased  by  simply  replacing  the  initial  mesh  with  a  mesh  obtained  after  one  subdi- 
vision step  applied  to  the  initial  mesh.  This  new  mesh  has  more  number  of  vertices 
but  defines  the  same  limit  surface.  Therefore,  the  dynamic  Catmull-Clark  surface 
model  can  be  subdivided  globally  to  increase  the  number  of  vertices  (control  points) 
of  the  model.  For  example,  after  one  step  of  global  subdivision,  the  initial  degrees  of 
freedom  p  (refer  to  Eqn.3.15  and  Eqn.3.16)  in  the  dynamic  system  will  be  replaced 
by  a  larger  number  of  degrees  of  freedom  q,  where  q  =  Ap.  A  is  a  global  subdivi- 
sion matrix  of  size  (M,N)  whose  entries  are  uniquely  determined  by  Catmull-Clark 
subdivision  rules  (see  Section  3.1  for  the  details  about  the  rules).  Thus,  p,  expressed 
as  a  function  of  q,  can  be  written  as 

p  =  (ArA)_1Arq  =  Afq,  (3.24) 
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where  A*  =  (ArA)~  AT.  Therefore,  Eqn.3.15  and  Eqn.3.16  can  be  rewritten  as 

s  =  (JAf)q,  (3.25) 

and 

s(ti,u,q)  =  (JAt)q,  (3.26) 

respectively.  Now  the  equation  of  motion  needs  to  be  derived  for  this  new  subdivided 
model  involving  a  larger  number  of  control  vertices  namely  q.  The  mass,  damping 
and  stiffness  matrices  need  to  be  recomputed  for  this  "finer"  level.  The  structure  of 
the  motion  equation  as  given  by  Eqn.2.5  remains  unchanged,  but  the  dimensionality 
and  the  entries  of  M,  D,K,p  and  fp  change  correspondingly  in  this  newly  obtained 
subdivided  level.  In  particular  the  motion  equation,  explicitly  expressed  as  a  function 
of  q,  can  be  written  as 

M,q  +  D,q  +  K,q  =  f„  (3.27) 

where  M,  =  J  J  fi(A^)  3T3A^dudv  and  the  derivation  of  D„  K,  and  f,  follow  suit. 
It  may  be  noted  that  further  subdivision,  if  necessary,  can  be  carried  out  in  a 
similar  fashion.  Therefore,  multilevel  dynamics  is  achieved  through  recursive  subdi- 
vision on  the  initial  set  of  control  vertices.  Users  can  interactively  choose  the  level  of 
detail  representation  of  the  dynamic  model  as  appropriate  for  their  modeling  and  de- 
sign requirements.  Alternatively,  the  system  can  automatically  determine  the  level  of 
subdivision  most  suitable  for  an  application  depending  on  some  application-specific 
criteria. 
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3.3     Finite  Element  Implementation 

The  evolution  of  the  generalized  coordinates  for  our  new  dynamic  surface  model 
can  be  determined  by  the  second-order  differential  equation  as  given  by  Eqn.2.5.  An 
analytical  solution  of  the  governing  differential  equation  can  not  be  derived  in  general. 
However,  an  efficient  numerical  implementation  can  be  obtained  using  the  finite  ele- 
ment method  [42].  The  limit  surface  of  the  dynamic  Catmull-Clark  subdivision  model 
is  a  collection  of  bicubic  patch  elements.  We  use  two  types  of  finite  elements  for  this 
purpose  -  normal  elements  (bicubic  patches  assigned  to  the  normal  faces  of  the  initial 
mesh)  and  special  elements  (collection  of  infinite  number  of  bicubic  patches  assigned 
to  each  extraordinary  vertex  of  the  initial  mesh).  The  shape  functions  for  the  nor- 
mal element  are  the  basis  functions  of  the  bicubic  surface  patch,  whereas  the  shape 
functions  for  the  special  element  are  the  collection  of  basis  functions  corresponding 
to  the  bicubic  patches  in  the  special  element.  In  the  current  implementation,  the 
M,  D  and  K  matrices  for  each  individual  normal  and  special  elements  are  calculated 
and  they  can  be  assembled  into  the  global  M,  D  and  K  matrices  that  appear  in  the 
corresponding  discrete  equation  of  motion.  In  practice,  we  never  assemble  the  global 
matrices  explicitly  in  the  interest  of  time  performance.  The  detailed  implementation 
is  explained  in  the  following  sections. 

3.3.1     Data  Structures 

The  limit  surface  of  the  dynamic  Catmull-Clark  subdivision  model  is  a  collection 
of  bicubic  patches,  and  hence  requires  us  to  keep  track  of  the  control  polygons  defining 
such  patches.  We  can  get  the  control  polygons  for  the  normal  elements  in  the  limit 
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-  list  of  vertices 

-  list  of  normal  elements 


Figure  3.6.  The  data  structure  used  for  dynamic  subdivision  surface  implementation. 


surface  from  the  initial  control  mesh  itself.  However,  we  need  to  locally  subdivide  the 
initial  control  mesh  around  the  extraordinary  vertices  to  obtain  the  control  polygons 
of  the  bicubic  patches  in  the  special  elements  on  the  limit  surface. 

A  subdivision  surface  defined  by  a  control  mesh  at  any  level  is  designed  as  a 
class  which  has  a  pointer  to  its  parent  mesh,  a  set  of  pointers  to  its  offspring  meshes 
(arising  out  of  local  subdivision  around  the  extraordinary  vertices  at  that  level),  a  list 
of  faces,  edges,  vertices  and  normal  elements  .  Face,  edge,  vertex  and  normal  elements 
are,  in  turn,  classes  which  store  all  the  connectivity  and  other  information  needed  to 
either  enumerate  all  the  patches  or  locally  subdivide  around  an  extraordinary  vertex 
in  that  level.  The  implementation  takes  the  initial  mesh  as  the  base  subdivision 
surface  object  (with  its  parent  pointer  set  to  NULL)  and  locally  subdivides  the  initial 
mesh  upto  a  user-defined  maximum  level  around  each  extraordinary  vertex  to  create 
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offspring  objects  at  different  levels  (Fig.3.6).  At  this  point,  let's  take  a  closer  look  at 
the  normal  and  special  element  data  structures  and  computation  of  the  corresponding 
local  M,  D  and  K  matrices. 

3.3.2  Normal  Elements 

Each  normal  element  is  a  bicubic  surface  patch  and  hence  is  defined  by  16 
vertices  (from  the  8-connected  neighborhood  of  the  corresponding  normal  face)  in  the 
initial  control  mesh.  Each  normal  element  keeps  a  set  of  pointers  to  those  vertices 
of  the  initial  mesh  which  act  as  control  points  for  the  given  element.  For  a  normal 
element,  the  mass,  damping  and  stiffness  matrices  are  of  size  (16, 16)  and  can  be 
computed  exactly  by  carrying  out  the  necessary  integrations  analytically.  The  matrix 
J  in  Eqn.3.18,  3.20  and  3.22  needs  to  be  replaced  by  Jp  (of  Eqn.3.2)  for  computation 
of  the  local  M,  D  and  K  matrices  respectively  of  the  corresponding  normal  element. 

3.3.3  Special  Elements 

Each  special  element  consists  of  an  infinite  number  of  bicubic  patches  in  the 
limit.  We  have  already  described  a  recursive  enumeration  of  the  bicubic  patches  of 
a  special  element  in  Section  3.2.2.  Let  us  now  consider  an  arbitrary  bicubic  patch  of 
the  special  element  in  some  level  j.  The  mass  matrix  Ms  of  this  patch  can  be  written 
as 

M,  =  njMpOs  (3.28) 

where  Mp  is  the  normal  element  mass  matrix  (scaled  by  a  factor  of  -^  to  take  into 
account  of  the  area  shrinkage  in  bicubic  patches  at  higher  level  of  subdivision)  and 
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ft,  is  the  transformation  matrix  of  the  control  points  of  that  arbitrary  patch  from  the 
corresponding  control  points  in  the  initial  mesh.  The  damping  and  stiffness  matrices 
for  the  given  bicubic  patch  can  be  derived  in  a  similar  fashion.  These  mass,  damping 
and  stiffness  matrices  from  various  levels  of  (local)  subdivision  can  then  be  assembled 
to  form  the  mass,  damping  and  stiffness  matrices  of  the  special  element.  It  may  be 
noted  that  the  stiffness  energy  due  to  plate  terms  diverges  at  the  extraordinary  points 
on  the  limit  surface.  We  solve  the  problem  using  a  slightly  different  approach  than 
the  one  used  in  Halstead  et  al.  [38].  When  the  area  of  the  bicubic  patch  obtained  via 
local  subdivision  of  the  initial  mesh  around  an  extraordinary  vertex  becomes  smaller 
than  the  display  resolution,  the  contribution  from  such  a  bicubic  patch  is  ignored  in 
computing  the  physical  matrices  of  the  corresponding  special  element.  The  number 
of  extraordinary  points  in  the  limit  surface  is  very  few,  and  the  above  mentioned 
approximation  is  found  to  work  well  in  practice. 

3.3.4     Force  Application 

The  force  f(u,  v,  t)  in  Eqn.3.23  represents  the  net  effect  of  all  externally  applied 
forces.  The  current  implementation  supports  spring,  inflation  as  well  as  image-based 
forces.  However,  other  types  of  forces  like  repulsion  forces,  gravitational  forces  etc. 
can  easily  be  implemented. 

To  apply  spring  forces,  a  spring  of  stiffness  k  can  be  connected  from  a  point  do 
to  a  point  x0  on  the  limit  surface  (or,  to  the  j-th  level  approximation  mesh),  the  net 
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applied  spring  force  being 

f  (x,  t)=  f       k(d0  -  s(x,  t))S(x  -  x0)dx,  (3.29) 

Jxesi 

where  8  is  the  unit  impulse  function  implying  f(x0,  t)  =  fc|do  —  s(x0,  t)\  and  vanishes 
elsewhere  on  the  surface.  However,  the  5  function  can  be  replaced  with  a  smooth 
kernel  to  spread  the  force  over  a  greater  portion  on  the  surface.  The  spring  forces  can 
be  applied  interactively  using  the  computer  mouse  or  the  points  from  which  forces 
need  to  be  applied  can  be  read  in  from  a  file. 

To  recover  shapes  from  3D  image  data,  we  synthesize  image-based  forces.  A  3D 
edge  detection  is  performed  on  a  volume  data  set  using  the  3D  Monga-Deriche(MD) 
operator  [69]  to  produce  a  3D  potential  field  P(x,y,z),  which  we  use  as  an  external 
potential  for  the  model.  The  force  distribution  is  then  computed  as 


f<*-y-*)=Aiivp(*,Mir  (3'30) 


where  A  controls  the  strength  of  the  force.  The  applied  force  on  each  vertex  at  the  j- 
th  approximation  level  is  computed  by  trilinear  interpolation  for  evaluating  Eqn.3.23 
in  Cartesian  coordinates.  More  sophisticated  image-based  forces  which  incorporate 
region-based  information  such  as  gradients  of  a  thresholded  fuzzy  voxel  classification 
can  also  be  used  to  yield  better  and  more  accurate  shape  recovery.  It  may  be  noted 
that  we  can  apply  spring  forces  in  addition  with  the  image-based  forces  by  placing 


54 

points  on  the  boundary  of  the  region  of  interest  in  each  slices  of  the  3D  volume  (MR, 

CT  etc.)  image  data. 

3.3.5     Discrete  Dynamic  Equation 

The  differential  equation  given  by  Eqn.2.5  is  integrated  through  time  by  dis- 
cretizing  the  time  derivative  of  p  over  time  steps  At.  The  state  of  the  dynamic 
subdivision  surface  at  time  t  +  At  is  integrated  using  prior  states  at  time  t  and 
t  -  At.  An  implicit  time  integration  method  is  used  in  the  current  implementation 
where  discrete  derivatives  of  p  are  calculated  using 


W  +  ±t)**t  +  At)~2PV  +  Pit-Atl  (3.3D 


and 

*ftj.AA      P(t  +  M)-p(t-At) 

p(f  +  At)  = ^ .  (3.32) 

The  elemental  mass,  damping  and  stiffness  matrices  can  be  assembled  to  get  the 
global  mass,  damping  and  stiffness  matrix  for  the  smooth  subdivision  surface  model. 
However,  we  do  not  assemble  these  global  sparse  matrices  explicitly  for  efficiency 
reasons.  For  the  time  varying  stiffness  matrix,  we  recompute  K  at  each  time  step. 
Using  Eqn.2.5,  Eqn.3.31  and  Eqn.3.32,  the  discrete  equation  of  motion  is  obtained  as 

(2M  +  DAt  +  2At2K)p(t  +  At)  =  2At%(t  +  At)  +  (D At  -  2M)p(t  -  At)  +  4Mp(t) . 

(3.33) 
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Figure  3.7.   Model  subdivision  to  increase  the  degrees  of  freedom  :  (a)  evolution  of 
the  initial  mesh,  and  (b)  the  corresponding  limit  surface  evolution  perceived  by  the 


This  linear  system  of  equations  is  solved  iteratively  between  each  time  step  using  the 
conjugate  gradient  method  [34,  75]. 

3.3.6     Model  Subdivision 

The  initialized  model  grows  dynamically  according  to  the  equation  of  motion 
(Eqn.2.5).  The  degrees  of  freedom  of  the  initialized  model  is  equal  to  the  number 
of  control  vertices  in  the  initial  mesh  as  mentioned  earlier.  When  an  equilibrium  is 
achieved  for  the  model,  the  number  of  control  vertices  can  be  increased  by  replacing 
the  original  initial  mesh  by  a  new  initial  mesh  obtained  through  one  step  of  Catmull- 
Clark  subdivision.  This  increases  the  number  of  degrees  of  freedom  to  represent  the 
same  (deformed)  smooth  limit  surface  and  a  new  equilibrium  position  for  the  model 
can  be  obtained.  This  process  is  depicted  schematically  in  Fig. 3. 7.  Model  subdivision 
might  be  needed  to  obtain  a  very  localized  effect  on  a  smooth  limit  surface.    For  a 
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shape  recovery  application,  one  may  start  with  a  very  simple  initial  model,  and  when 
an  approximate  shape  is  recovered,  the  degrees  of  freedom  can  be  increased  to  obtain 
a  new  equilibrium  position  for  the  model  with  a  better  fit  to  the  given  data  set.  The 
error  of  fit  criteria  for  the  discrete  data  is  based  on  distance  between  the  data  points 
and  the  points  on  the  limit  surface  where  the  corresponding  springs  are  attached. 
In  the  context  of  image-based  forces,  if  the  model  energy  does  not  change  between 
successive  iterations  indicating  an  equilibrium  for  the  given  resolution,  the  degrees  of 
freedom  for  the  model  can  be  increased  by  the  above-mentioned  replacement  scheme 
until  the  model  energy  is  sufficiently  small  and  the  change  in  model  energy  between 
successive  iterations  becomes  less  than  a  pre-specified  tolerance. 

3.4     Generalization  of  the  Approach 

The  proposed  approach  can  be  generalized  for  other  approximating  subdivision 
schemes.  However,  a  more  general  approach  is  presented  in  Chapter  5  for  deriving  the 
dynamic  framework.  Dynamic  Catmull-Clark  subdivision  surface  model  is  reformu- 
lated in  Section  5.2  using  this  general  approach.  A  dynamic  framework  for  another 
popular  approximating  subdivision  scheme  namely,  Loop's  subdivision  scheme,  is  also 
discussed  in  Chapter  5. 


CHAPTER  4 
DYNAMIC  BUTTERFLY  SUBDIVISION  SURFACES 


In  the  previous  chapter,  a  dynamic  framework  has  been  presented  for  an  approx- 
imating subdivision  scheme  namely,  Catmull-Clark  subdivision  scheme.  In  this  chap- 
ter, a  dynamic  framework  is  developed  for  (modified)  butterfly  subdivision  scheme, 
which  is  an  interpolatory  subdivision  technique.  First,  a  brief  overview  of  the  (mod- 
ified) butterfly  subdivision  scheme  is  presented.  Next,  a  local  parameterization  tech- 
nique for  the  limit  surface  obtained  via  (modified)  butterfly  subdivision  is  discussed. 
This  parameterization  scheme  is  then  used  to  derive  the  dynamic  model.  The  imple- 
mentation details  are  also  discussed. 

4.1     Overview  of  the  (Modified)  Butterfly  Subdivision  Scheme 

The  butterfly  subdivision  scheme  [30],  like  many  other  subdivision  schemes  used 
in  geometric  design  literature/applications,  starts  with  an  initial  triangular  mesh 
which  is  also  known  as  the  control  mesh.  The  vertices  of  the  control  mesh  are  known 
as  the  control  points.  In  each  step  of  subdivision,  the  initial  (control)  mesh  is  refined 
through  the  transformation  of  each  triangular  face  into  a  patch  with  four  smaller 
triangular  faces.  After  one  step  of  refinement,  the  new  mesh  in  the  finer  level  retains 
the  vertices  of  each  triangular  face  in  the  previous  level  and  hence,  interpolates  the 
coarser  mesh  in  the  previous  level.   In  addition,  every  edge  in  each  triangular  face 
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(a) 


(b) 


Figure  4.1.   (a)  The  control  polygon  with  triangular  faces;  (b)  mesh  obtained  after 
one  subdivision  step  using  butterfly  subdivision  rules. 


-W 


(a) 
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Figure  4.2.  (a)  The  contributing  vertices  in  the  j-th  level  for  the  vertex  in  the  j+1- 
th  level  corresponding  to  the  edge  between  v{  and  v32;  (b)  the  weighing  factors  for 
different  vertices. 
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is  spilt  by  adding  a  new  vertex  whose  position  is  obtained  by  an  affine  combination 
of  the  neighboring  vertex  positions  in  the  coarser  level.  For  instance,  the  mesh  in 
Fig.4.1(b)  is  obtained  by  subdividing  the  initial  mesh  shown  in  Fig.4.1(a)  once.  It 
may  be  noted  that  all  the  newly  introduced  vertices  corresponding  to  the  edges  in  the 
original  mesh  have  degree  6,  whereas  the  position  and  degree  of  the  original  vertices 
do  not  change  in  the  subdivided  mesh. 

In  the  original  butterfly  scheme,  the  new  vertices  corresponding  to  the  edges  in 
the  previous  level  are  obtained  using  an  eight-point  stencil  as  shown  in  Fig.4.2(a).  To 
show  how  this  stencil  is  used,  a  vertex  to  be  introduced  and  the  contributing  vertices 
from  the  stencil  are  highlighted  in  Fig.4.1(a).  The  name  of  the  scheme  originated 
from  the  "butterfly"-like  configuration  of  the  contributing  vertices.  The  weighing 
factors  for  different  contributing  vertex  positions  are  shown  in  Fig.4.2(b).  The  vertex 
ejj '  in  the  j  +  1-th  level  of  subdivision,  corresponding  to  the  edge  connecting  vertices 
Vj  and  v?  at  level  j,  is  obtained  by 

eft1  =  0.5(v{  +  v{)  +  2w{v{  +  vJ4)  -  w{vi  +  v£  +  vJ7  +  vg),  (4.1) 

where  0  <  W  <  1,  and  v'  denotes  the  position  of  the  i-th  vertex  at  the  j-th  level. 

The  butterfly  subdivision  scheme  produces  a  smooth  C1  surface  in  the  limit 
except  at  the  extraordinary  points  corresponding  to  the  extraordinary  vertices  (ver- 
tices with  degree  not  equal  to  6)  in  the  initial  mesh  [109].  All  the  vertices  introduced 
through  subdivision  have  degree  6,  and  therefore,  the  number  of  extraordinary  points 
in  the  smooth  limit  surface  equals  the  number  of  extraordinary  vertices  in  the  initial 
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Figure  4.3.  (a)  The  weighing  factors  of  contributing  vertex  positions  for  an  edge 
connecting  two  vertices  of  degree  6;  (b)  the  corresponding  case  when  one  vertex  is  of 
degree  n  and  the  other  is  of  degree  6. 
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mesh.  Recently,  this  scheme  has  been  modified  by  Zorin  et  al.  [109]  to  obtain  better 
smoothness  properties  at  the  extraordinary  points.  This  modified  scheme  is  known 
as  modified  butterfly  subdivision  technique.  In  Zorin  et  al.  [109],  all  the  edges  have 
been  categorized  into  three  classes  :  (i)  edges  connecting  two  vertices  of  degree  6  (a 
10  point  stencil,  as  shown  in  Fig.4.3(a),  is  used  to  obtain  the  new  vertex  positions 
corresponding  to  these  edges),  (ii)  edges  connecting  a  vertex  of  degree  6  and  a  vertex 
of  degree  n  ^  6  (the  corresponding  stencil  to  obtain  new  vertex  position  is  shown  in 
Fig.4.3(b),  where  q  =  .75  is  the  weight  associated  with  the  vertex  of  degree  n  ^  6, 
and  st  =  (0.25  +  cos(2iri/n)  +  O.bcos(4ni/n))/n,  i  =  0, 1, . . . ,  n  —  1,  are  the  weights 
associated  with  the  vertices  of  degree  6),  and  (iii)  edges  connecting  two  vertices  of 
degree  n  ^  6.  The  last  case  can  not  occur  except  in  the  initial  mesh  as  the  newly 
introduced  vertices  are  of  degree  6,  and  the  new  vertex  position  in  this  last  case  is 
obtained  by  averaging  the  positions  obtained  through  the  use  of  stencil  shown  in 
Fig.4.3(b)  at  each  of  those  two  extraordinary  vertices. 

4.2     Formulation 

In  this  section,  a  systematic  formulation  of  the  dynamic  framework  for  the  modi- 
fied butterfly  subdivision  scheme  is  presented.  Unlike  the  approximating  schemes,  the 
limit  surface  obtained  via  modified  butterfly  subdivision  does  not  have  any  analytic 
expression  even  in  a  regular  setting.  Therefore,  derivation  of  a  local  parameterization 
suitable  for  developing  the  dynamic  framework  for  the  limit  surface  is  the  key  issue 
here  which  is  discussed  next. 
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Figure  4.4.  The  smoothing  effect  of  the  subdivision  process  on  the  triangles  of  the 
initial  mesh. 


4.2.1     Local  Parameterization 

The  smooth  limit  surface  defined  by  the  modified  butterfly  subdivision  technique 
is  of  arbitrary  topology  where  a  global  parameterization  is  impossible.  Nevertheless, 
the  limit  surface  can  be  locally  parameterized  over  the  domain  defined  by  the  initial 
mesh  following  a  similar  approach  described  in  Lounsbery  et  al.  [53].  The  idea 
is  to  track  any  arbitrary  point  on  the  initial  mesh  across  the  meshes  obtained  via 
the  subdivision  process  (see  Fig.4.4  and  Fig.4.5),  so  that  a  correspondence  can  be 
established  between  the  point  being  tracked  in  the  initial  mesh  and  its  image  on  the 
limit  surface. 

The  modified  butterfly  subdivision  scheme  starts  with  an  initial  mesh  consisting 
of  a  set  of  triangular  faces.  The  recursive  application  of  the  subdivision  rules  smoothes 
out  each  triangular  face,  and  in  the  limit,  a  smooth  surface  is  obtained  which  can  also 
be  considered  as  a  collection  of  smooth  triangular  patches.  The  subdivision  process 
and  the  triangular  decomposition  of  the  limit  surface  is  depicted  in  Fig.4.4.  Note 
that,  the  limit  surface  can  be  represented  by  the  same  number  of  smooth  triangular 
patches  as  that  of  the  triangular  faces  in  the  initial  mesh.  Therefore,  the  limit  surface 
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s  can  be  expressed  as 


s  =  E^.  (4-2) 


where  n  is  the  number  of  triangular  faces  in  the  initial  mesh  and  St  is  the  smooth 
triangular  patch  in  the  limit  surface  corresponding  to  the  fc-th  triangular  face  in  the 
initial  mesh. 

The  stage  is  now  set  for  describing  the  parameterization  of  the  limit  surface  over 
the  initial  mesh.  The  process  is  best  explained  through  the  following  example.  A 
simple  planar  mesh  shown  in  Fig.4.5(a)  is  chosen  as  the  initial  mesh.  An  arbitrary 
point  x  inside  the  triangular  face  abc  is  tracked  over  the  meshes  obtained  through 
subdivision.  The  vertices  in  the  initial  mesh  are  darkly  shaded  in  Fig.4.5.  After 
one  step  of  subdivision,  the  initial  mesh  is  refined  by  addition  of  new  vertices  which 
are  lightly  shaded.  Another  subdivision  step  on  this  refined  mesh  leads  to  a  finer 
mesh  with  introduction  of  new  vertices  which  are  unshaded.  It  may  be  noted  that 
any  point  inside  the  smooth  triangular  patch  in  the  limit  surface  corresponding  to  the 
face  abc  in  the  initial  mesh  depends  only  on  the  vertices  in  the  initial  mesh  which 
are  within  the  2-neighborhood  of  the  vertices  a,  b  and  c  due  to  the  local  nature  of 
the  subdivision  process.  For  example,  the  vertex  d,  introduced  after  first  subdivision 
step,  can  be  obtained  using  the  10  point  stencil  shown  in  Fig.4.3(a)  on  the  edge  ab. 
All  the  contributing  vertices  in  the  initial  mesh  are  within  the  1-neighborhood  of  the 
vertices  a  and  b.  A  10  point  stencil  can  be  used  again  in  the  next  subdivision  step 
on  the  edge  db  to  obtain  the  vertex  g.  Some  of  the  contributing  vertices  at  this  level 
of  subdivision,  for  example,  the  (lightly  shaded)  1-neighbors  of  the  vertex  b  (except 
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Figure  4.5.  Tracking  a  point  x  through  various  levels  of  subdivision  :  (a)  initial 
mesh,  (b)  the  selected  section  (enclosed  by  dotted  lines)  of  the  mesh  in  (a),  after  one 
subdivision  step,  (c)  the  selected  section  of  the  mesh  in  (b),  after  another  subdivision 
step. 
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d  and  e)  in  Fig.4.5(b),  depend  on  some  vertices  in  the  initial  mesh  which  are  within 
the  2-neighborhood  of  the  vertices  a,  b  and  c  in  the  initial  mesh. 

In  the  rest  of  the  formulation,  superscripts  are  used  to  indicate  the  subdivision 
level.  For  example,  v}uvw  denotes  the  collection  of  vertices  at  level  j  which  control 
the  smooth  patch  in  the  limit  surface  corresponding  to  the  triangular  face  uvw  at 
the  j'-th  level  of  subdivision.  Let  vjj^  be  the  collection  of  vertices  in  the  initial  mesh 
which  are  within  the  2-neighborhood  of  the  vertices  a,  b  and  c  (marked  black  in 
Fig.4.5(a)).  Let  the  number  of  such  vertices  be  r.  Then,  the  vector  v^,.,  which  is  the 
concatenation  of  the  (x,  y,  z)  positions  for  all  the  r  vertices,  is  of  dimension  3r.  These 
r  vertices  control  the  smooth  triangular  patch  in  the  limit  surface  corresponding  to 
the  triangular  face  abe  in  the  initial  mesh.  Now,  there  exists  four  subdivision  matrices 
(Aotc)(,  (Ao6c)(,  (Aabc)r  and  (Aate)ro  of  dimension  (3r,3r)  such  that 

vldf     =     (Aa6c),V°6c, 

VL   =    (A«te),vJ^., 

vi,,    =    (Ao4c)mvl,  (4.3) 

where  the  subscripts  t,  I,  r  and  m  denote  top,  left,  right  and  middle  triangle  positions 
respectively  (indicating  the  relative  position  of  the  new  triangle  with  respect  to  the 
original  triangle),  and  v^,  vjed,  vL,  and  v^  are  the  concatenation  of  the  (x,  y,  z) 
positions  for  the  vertices  in  the  2-neighborhood  of  the  corresponding  triangle  in  the 
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newly  obtained  subdivided  mesh.  The  new  vertices  in  this  level  of  subdivision  are 
lightly  shaded  in  Fig.4.5(b).  The  2-neighborhood  configuration  of  the  vertices  in  the 
newly  obtained  triangles  is  exactly  the  same  as  that  of  the  original  triangle,  hence 
local  subdivision  matrices  are  square  and  the  vector  dimensions  on  both  sides  of 
Eqn.4.3  are  the  same.  This  concept  is  further  illustrated  in  Fig.4.6. 

Carrying  out  one  more  level  of  subdivision,  a  new  set  of  vertices  which  are 
unshaded  in  Fig.4.5(c)  are  obtained  along  with  the  old  vertices.  Adopting  a  similar 
approach  as  in  the  derivation  of  Eqn.4.3,  it  can  be  shown  that 


vl, 

= 

(Abed 

tVl.d 

vlhg 

= 

(A-bed 

l^led 

vlh 

= 

(Ab^ 

rVled 

VJW 

= 

(Afed 

mVbed 

(4.4) 


The  relative  position  of  the  triangular  face  dgi  in  Fig.4.5(c)  with  respect  to  the 
triangular  face  bed  is  topologically  the  same  as  of  the  triangular  face  adf  in  Fig.4.5(b) 
with  respect  to  the  triangular  face  abc.  Therefore,  one  can  write  (A^),  =  {Aabc)r 
Using  similar  reasoning,  Eqn.4.4  can  be  rewritten  as 

vdgi      =      {Abed)tvled  =  {A.abc)tvled 
Vlhg      =      (A(,ed),vL  =  (Aabc)ybed 

VL    =    {Abed)rvlbed  =  (A0i,c)rv»ed 
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(a) 


(l.) 


(c) 


(d) 


Figure  4.6.  Different  set  of  control  points  at  a  subdivided  level  obtained  by  applying 
different  subdivision  matrices  on  a  given  set  of  control  points  in  a  coarser  mesh. 
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vSfc.      =      (A*erf)mvLd  =  (Aofc),^,,.  (4.5) 

Combining  Eqn.4.3  and  Eqn.4.5,  it  can  be  shown  that 
v29,    =    (A„6c)((Aoi,c),v26c, 

VL    =    (A06C)r(Aajc),v2jc, 

vs\,    =    (Aa6c)m(Aatc)1vl.  (4.6) 

Let  x  be  a  point  with  barycentric  coordinates  (a°bc,  0^,  7°(,c)  inside  the  trian- 
gular face  abc.  When  the  initial  mesh  is  subdivided,  x  becomes  a  point  inside  the 
triangular  face  bed  with  barycentric  coordinates  (a^,  0^,  7jerf).  Another  level  of 
subdivision  causes  x  to  be  included  in  the  triangular  face  dgi  with  barycentric  coor- 
dinates (o?dgi,0dsi>ldgi)-  Let  s3abc  denote  the  j-th  level  approximation  of  the  smooth 
triangular  patch  satc  in  the  limit  surface  corresponding  to  the  triangular  face  abc  in 
the  initial  mesh.  Now  vj^  can  be  written  as 


v2tc  =  iax,  bx,  cx, . . .,  ay,  by,  Cy, . . .,  a„  bz,  cz, . .  )T 

where  the  subscripts  x,  y  and  z  indicate  the  x,  y  and  z  coordinates  respectively  of  the 
corresponding  vertex  position.  The  expressions  for  Vjerf  and  \dgl  can  also  be  written 


70 


in  a  similar  manner.  Next,  the  matrix  B°tc  can  be  constructed  as  follows: 


Bl(x) 


0^, /&»-&*.  <>■•■•, 0,0,. ..,0,0,---0 


0,...,0,al,/3l,7L,0,...,0,0,...,0 


0,...,0,0,...,0,al,/?l,7l,0,...,0 


The  matrices  Bjeii  and  B^9i  can  also  be  constructed  in  a  similar  fashion.  Now  s°,,c(x), 
sjj^x),  and  s^c(x)  can  be  written  as 


cw  ■ 

-    B°tc(x)vl, 

4tcW   ■ 

=  BL(x)vL    - 

=    Bt^xJfA^JjvJh 

»LM  ■ 

=    B5si(x)v^ 

=     B^,(x)(Aai,c)(vtled 

(4.7) 


Proceeding  in  a  similar  way,  the  expression  for  s^c(x),  j-th  level  approximation 
of  sa6c(x),  is  given  by 


■LW 


Bi„w(x)(Aotc)m...(Aa6c)((Ao6c),vl 
B„„ro(x)(Aa6(.)vo(,c 


(4.8) 


where  x  is  inside  the  triangular  face  uvw  at  level  j  (with  an  assumption  that  uvw  is  the 
triangular  face  in  the  middle  with  respect  to  its  coarser  level  original  triangular  face  in 
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the  previous  level),  (AJJ  =  (Aabc)m  ■  ■  ■  {Aabc)t(Aabc)l  and  B^(x)  =  Bjl)l„(x)(A^c). 
It  may  be  noted  that  the  sequence  of  applying  (Aabc)t,  (Ao6c)„  (Aabc)r  and  (Aabc)m 
depends  on  the  triangle  inside  which  the  tracked  point  x  falls  after  each  subdivision 
step.  Finally,  the  local  parameterization  process  can  be  completed  by  writing 


satc(x)  =  ( Jim  BJQiK(x))vl  =  Bo(,c(x)vl.  (4.9) 

J-fOO 


In  the  above  equation,  Ba(,c  is  the  collection  of  basis  functions  at  the  vertices  of 
v®bc.  It  may  also  be  noted  that  the  modified  butterfly  subdivision  scheme  is  a  sta- 
tionary subdivision  process,  and  hence  new  vertex  positions  are  obtained  by  affine 
combinations  of  nearby  vertices.  This  guarantees  that  each  row  of  the  matrices 
{Aabc)ti  {Aabc)t,  {Aabc)T  and  (Aaj,,.)m  sums  to  one.  The  largest  eigenvalue  of  such  ma- 
trices is  1  and  therefore  the  limit  in  Eqn.4.9  exists.  Now,  assuming  the  triangular 
face  abc  is  the  fc-th  face  in  the  initial  mesh,  Eqn.4.9  can  be  rewritten  as 

S*(x)  =  B*(x)v°  =  B*(x)A4p,  (4.10) 

where  p  is  the  concatenation  of  the  (x,y,z)  positions  of  all  the  vertices  in  the  initial 
mesh  and  the  matrix  Ak,  when  post-multiplied  by  p,  selects  the  vertices  v°  controlling 
the  k-th  smooth  triangular  patch  in  the  limit  surface.  If  there  are  t  vertices  in  the 
initial  mesh  and  r  of  them  control  the  fc-th  patch,  then  p  is  a  vector  of  dimension  3£, 
At  is  a  matrix  of  dimension  (3r,  3f),  and  B/i(x)  is  a  matrix  of  dimension  (3,3r). 
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Combining  Eqn.4.2  and  Eqn.4.10,  it  can  be  shown  that 


s(x)  =  (£  Bfe(x)A,)p  =  J(x)p,  (4.11) 


where  J,  a  matrix  of  dimension  (3,3t),  is  the  collection  of  basis  functions  for  the 
corresponding  vertices  in  the  initial  mesh.  The  vector  p  is  also  known  as  the  degrees 
of  freedom  vector  of  the  smooth  limit  surface  s. 

4.2.2     Dynamics 

An  expression  of  the  limit  surface  obtained  via  modified  butterfly  subdivision 
process  is  derived  in  the  last  section.  Now  a  dynamic  framework  for  the  limit  surface 
can  be  derived  using  this  expression.  The  derivation  of  the  dynamic  model  from  this 
point  is  exactly  similar  in  spirit  as  that  of  the  dynamic  Catmull-Clark  subdivision 
model  presented  in  Section  3.2.4.  The  vertex  positions  in  the  initial  mesh  defining 
the  smooth  limit  surface  s  are  treated  as  a  function  of  time  in  order  to  embed  the 
modified  butterfly  subdivision  model  in  a  dynamic  framework.  The  velocity  of  this 
surface  model  can  be  expressed  as 

s(x,p)=J(x)p,  (4.12) 

where  an  overstruck  dot  denotes  a  time  derivative  and  x  6  5°,  S°  being  the  domain 
defined  by  the  initial  mesh.  As  in  the  case  of  the  dynamic  Catmull-Clark  subdivision 
surfaces,  the  Lagrangian  equation  of  motion  (Eqn.2.1)  is  used  to  derive  the  dynamic 
modified  butterfly  subdivision  surface  model. 


73 


Let  fi  be  the  mass  density  function  of  the  surface.  Then  the  kinetic  energy  of 
the  surface  is 

T    =    \f       /i(x)sT(x)s(x)dx      =       VMp,  (4.13) 

where  (using  Eqn.4.12)  M  =  /x6So  /j(x)JT(x)J(x)dx  is  the  mass  matrix  of  dimen- 
sion (3i,3i).  Similarly,  let  7  be  the  damping  density  function  of  the  surface.  The 
dissipation  energy  is 

F    =    \f       7(x)sT(x)s(x)dx      =       JpTDp,  (4.14) 

2  Jxeso  2 

where  D  =  JX€S0  7(x)JT(x)J(x)dx  is  the  damping  matrix  of  dimension  (3i,3t).  The 
potential  energy  of  the  smooth  limit  surface  can  be  expressed  as 


U    =    ipTKp,  (4.15) 


where  K  is  the  stiffness  matrix  of  dimension  (3t,  3t)  obtained  by  assigning  various 
internal  energies  to  a  discretized  approximation  of  the  limit  surface  and  is  detailed 
in  Section  4.3. 

Using  the  expressions  for  the  kinetic,  dissipation  and  potential  energy  in  Eqn.2.1, 
the  same  motion  equation  as  in  Eqn.2.5  can  be  obtained.  The  generalized  force  vector 
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fp  can  be  obtained  is  a  similar  fashion  described  in  Section  3.2.4  and  is  given  by 

f.=  f       Jr(x)f(x,t)dx.  (4.16) 

Jx€S° 

4.2.3     Multilevel  Dynamics 

The  initial  mesh  of  the  dynamic  modified  butterfly  subdivision  surface  model 
can  be  subdivided  to  increase  the  degrees  of  freedom  for  model  representation.  The 
situation  is  exactly  similar  to  that  of  the  dynamic  Catmull-Clark  subdivision  surface 
model  as  described  in  Section  3.2.5  and  a  similar  example  is  used  for  illustration. 
After  one  step  of  modified  butterfly  subdivision,  the  initial  degrees  of  freedom  p  (refer 
to  Eqn.4.11  and  Eqn.4.12)  in  the  dynamic  system  is  replaced  by  a  larger  number  of 
degrees  of  freedom  q,  where  q  =  Bp.  B  is  a  global  subdivision  matrix  of  size  (3s,  3t) 
whose  entries  are  uniquely  determined  by  the  weights  used  in  the  modified  butterfly 
subdivision  scheme  (see  Section  4.1  for  the  weights).  Thus,  p,  expressed  as  a  function 
of  q,  can  be  written  as 

p  =  (BTB)~lBrq  =  Btq,  (4.17) 

where  B+  =  (BTBf1BT.  Therefore,  Eqn.4.11  and  Eqn.4.12  is  modified  as 

s(x)  =  (J(x)Bt)q,  (4.18) 

and 

s(x,q)  =  (J(x)Bt)q,  (4.19) 
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Figure  4.7.  (a)  An  initial  mesh,  and  (b)  the  corresponding  limit  surface.  The  domains 
of  the  shaded  elements  in  the  limit  surface  are  the  corresponding  triangular  faces  in 
the  initial  mesh.  The  encircled  vertices  in  (a)  are  the  degrees  of  freedom  for  the 
corresponding  element. 


respectively.    The  motion  equation,  explicitly  expressed  as  a  function  of  q,  can  be 
written  as 

M,q  +  D,q  +  K„q  =  f„  (4.20) 


where  M,  =  /x€Si  /i(x)(Bt)  JT(x)J(x)Btrfx,  S'1  being  the  domain  defined  by  the 
newly  obtained  subdivided  mesh.  The  derivations  of  D,,  K,  and  f,  are  similar. 
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4,3     Finite  Element  Implementation 

In  the  previous  section,  the  smooth  limit  surface  obtained  via  modified  butterfly 
subdivision  process  is  expressed  as  a  function  of  the  control  vertex  positions  in  its 
initial  mesh;  mass  and  damping  distribution,  internal  deformation  energy  and  forces 
are  assigned  to  the  limit  surface  in  order  to  develop  the  corresponding  physical  model. 
In  this  section,  the  implementation  of  this  physical  model  is  detailed.  It  may  be  noted 
that  even  though  the  formulation  of  the  dynamic  framework  for  modified  butterfly 
subdivision  scheme  is  similar  to  that  of  Catmull-Clark  subdivision  scheme  once  the 
local  parameterization  is  derived,  the  finite  element  implementation  is  totally  different 
as  the  limit  surface  does  not  have  any  analytic  expression  in  case  of  the  modified 
butterfly  subdivision  scheme. 

In  Section  4.2  it  was  pointed  out  that  the  smooth  limit  surface  obtained  by  the 
recursive  application  of  the  modified  butterfly  subdivision  rules  can  be  represented  by 
a  set  of  smooth  triangular  patches,  each  of  which  is  represented  by  a  finite  element. 
The  shape  (basis)  function  of  this  finite  element  is  obtained  by  smoothing  a  hat 
function  through  repeated  application  of  the  modified  butterfly  subdivision  rules. 
The  number  of  finite  elements  in  the  smooth  limit  surface  is  equal  to  the  number 
of  triangular  faces  in  the  initial  mesh  as  mentioned  earlier  (refer  Fig. 4. 4  and  4.7). 
A  detailed  discussion  on  how  to  derive  the  mass,  damping  and  stiffness  matrices 
for  these  elements  is  presented  next.  These  elemental  matrices  can  be  assembled 
to  obtain  the  global  physical  matrices  M,  D  and  K,  and  a  numerical  solution  to 
the  governing  second-order  differential  equation  as  given  by  Eqn.2.5  can  be  obtained 
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using  finite  element  analysis  techniques  [42].  The  same  example  as  in  Section  4.2 
(refer  Fig.4.5)  is  used  to  develop  the  related  concepts.  The  concept  of  elements  along 
with  the  control  vertices  and  their  corresponding  domains  is  further  illustrated  in 
Fig.4.7. 

Now  it  will  be  shown  how  to  derive  the  mass,  damping  and  stiffness  matrices 
for  the  element  corresponding  to  the  triangular  face  abc  in  Fig.4.5.  The  derivations 
hold  for  any  element  in  general. 

4.3.1     Elemental  Mass  and  Damping  matrices 

The  mass  matrix  for  the  element  given  by  s0(,c  (Eqn.4.9)  can  be  written  as 

Mo6c=/        /i(x){Botc(x)}r{Bo(K(x)}rfx.  (4.21) 

However,  from  Eqn.4.9  it  is  known  that  the  basis  functions  corresponding  to  the 
vertices  in  v°bc  which  are  stored  as  entries  in  B„;,c  are  obtained  as  a  limiting  process. 
These  basis  functions  do  not  have  any  analytic  form,  hence  computing  the  inner 
product  of  such  basis  functions  as  needed  in  Eqn.4.21  is  a  challenging  problem.  In 
Lounsbery  et  al.  [53],  an  outline  is  provided  on  the  computation  of  these  inner 
products  without  performing  any  integration.  In  this  dissertation,  a  different  yet 
simpler  approach  is  developed  to  solve  this  problem.  The  smooth  triangular  patch  in 
the  limit  surface  corresponding  to  the  face  abc  in  the  initial  mesh  is  approximated  by 
a  triangular  mesh  with  4J  faces  obtained  after  j  levels  of  subdivision  of  the  original 
triangular  face  abc  (each  subdivision  step  splits  one  triangular  face  into  4  triangular 
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faces).  Then  the  mass  matrix  can  be  expressed  as 

M„6c  =  W       Mx){BJatc(x)}r{B^c(x)}rfx.  (4.22) 

The  j-th  level  approximation  of  the  corresponding  basis  functions  can  be  explicitly 
evaluated  (refer  Eqn.4.8  for  an  expression  of  B^,c).  An  important  point  to  note  is 
that  Eqn.4.8  involves  several  matrix  multiplications  and  hence  can  be  very  expensive 
to  evaluate.  However,  the  matrix  (A^,)(=  (Aok)m  . . .  (Ao6c)((Aa6c),)  in  Eqn.4.8 
encodes  how  vertices  in  the  2-neighborhood  of  the  triangular  face  uvw  at  level  j  are 
related  to  the  vertices  in  the  2-neighborhood  of  the  triangular  face  abc  in  the  initial 
mesh.  In  the  implementation,  how  a  new  vertex  is  obtained  from  the  contributing 
vertices  in  its  immediate  predecessor  level  is  tracked.  The  information  stored  in  ( A^c) 
can  be  obtained  by  tracking  the  contributions  from  level  j  to  level  0  recursively.  Thus 
the  entries  of  ( Ajabc)  can  be  determined  without  forming  any  local  subdivision  matrices 
and  thereby  avoiding  subsequent  matrix  multiplications. 

By  choosing  a  sufficiently  high  value  of  j,  a  reasonably  good  approximation  of 
the  elemental  mass  matrices  is  achieved.  The  computations  involved  in  evaluating 
the  integrals  in  Eqn.4.22  are  eliminated  by  using  discrete  mass  density  function  which 
has  non-zero  values  only  at  the  vertex  positions  of  the  j-th  subdivision  level  mesh. 
Therefore,  the  approximation  of  the  mass  matrix  for  the  element  can  be  written  as 


Motc  =  £MvD{B'atcK)}T{B^(vD},  (4.23) 
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where  k  is  the  number  of  vertices  in  the  triangular  mesh  with  V  faces.  This  approxi- 
mation has  been  found  to  be  very  effective  and  efficient  for  implementation  purposes. 
The  elemental  damping  matrix  can  be  obtained  in  a  similar  fashion. 

4.3.2     Elemental  Stiffness  Matrix 

The  internal  energy  is  assigned  to  each  element  in  the  limit  surface,  thereby 
defining  the  internal  energy  of  the  smooth  subdivision  surface  model.  A  similar 
approach  as  in  the  derivation  of  the  elemental  mass  and  damping  matrix  is  taken 
and  the  internal  energy  is  assigned  to  a  j-th  level  approximation  of  the  element. 

Three  types  of  internal  energy  are  used  -  tension,  stiffness  and  spring  energy. 
For  the  examples  used  throughout  the  paper,  this  energy  at  the  j-th  level  of  approx- 
imation can  be  written  as 

Eabc  *  E^  =  (Eibc)t  +  (EUm  +  (E^U-  <4-24' 

where  the  subscripts  t,  st  and  sp  denote  tension,  stiffness  and  spring  respectively. 

The  expression  for  the  tension  energy,  which  is  essentially  equivalent  to  the  first 
order  strain  (membrane)  energy  [96],  is 


=    \h{vlbf  (K{bc)t{<bch  (4-25) 


where  kt  is  a  constant,  v/  and  \'m,  the  /-th  and  m-th  vertex  in  the  j-th  level  mesh, 
are  in  the  1-neighborhood  of  each  other,  Q  is  the  domain  defined  by  all  such  vertex 
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pairs,  and  v{bc  is  the  concatenation  of  the  (x,y,z)  positions  of  all  the  vertices  in  the 
j-th  subdivision  level  of  the  triangular  face  abc  in  the  initial  mesh. 

Similarly,  the  expression  for  stiffness  energy,  which  is  equivalent  to  the  second 
order  strain  (thin  plate)  energy  [96],  can  be  written  as 

*■       n 
=    jM*i.}Tff2taU*W.  (4-26) 

where  v/,  \'m  and  v'n  are  the  three  vertices  of  a  triangular  face.  The  summation  in- 
volves three  terms  corresponding  to  each  triangular  face,  and  is  over  all  the  triangular 
faces  in  the  mesh  at  the  j-th  level  of  subdivision. 

Finally,  a  spring  energy  term  is  added  which  is  given  by 


j  _     1  ^  (felmUlvf-VJI-4.).    ,  _      ■ 

=    ^{vL}T(K^)spK6c},  (4-27) 


where  (klm)  is  the  spring  constant,  Q  is  the  domain  as  in  Eqn.4.25  and  ltm  is  the 
natural  length  of  the  spring  connected  between  \j  and  vJm.  It  may  be  noted  that 
the  entries  in  (Kjabc)  depend  on  the  distance  between  the  connected  vertices  and 
hence,  (K3abc)s  ,  unlike  other  elemental  matrices,  is  a  function  of  time  which  needs  to 
be  recomputed  in  each  time  step. 
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Combining  the  expressions  for  tension,  stiffness  and  spring  energy,  it  can  be 
shown  that 


1,    i    ,T 

2 


vi.}  {*t(Kjjt  +  fc,t(Ki6c)s(  +  (KJatc)sp}K6c} 


=    i{(AL){vl}}T(Klbc)(A^c){vl} 

=    i{vl}r{(A^)T(K^)(A^)}{vl},  (4.28) 


where  (AJjJ  and  v°abc  are  same  as  in  Eqn.4.8.    Therefore,  the  expression  for  the 
elemental  stiffness  matrix  is  given  by 


K«6c  =  (Al6/(KL)(AJa6c).  (4.29) 


It  may  be  noted  that  the  matrix  multiplications  for  constructing  K0j,c  are  avoided  in 
the  implementation  by  following  the  same  technique  described  in  Section  4.3.1. 

4.3.3     Other  Implementation  Issues 

The  techniques  of  applying  synthesized  forces  on  the  smooth  limit  surface  ob- 
tained via  modified  butterfly  subdivision  process  are  similar  to  those  described  in 
Section  3.3.4  in  the  context  of  dynamic  Catmull-Clark  subdivision  surface  model. 
The  techniques  of  model  subdivision  and  discrete  dynamic  equation  derivation  are 
also  the  same. 
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4.4     Generalization  of  the  Approach 

The  proposed  approach  of  deriving  the  dynamic  framework  for  modified  butter- 
fly subdivision  scheme  is  generalized  for  other  interpolatory  subdivision  schemes  in 
the  next  chapter,  Chapter  5. 


CHAPTER  5 
UNIFIED  DYNAMIC  FRAMEWORK  FOR  SUBDIVISION  SURFACES 


Dynamic  framework  for  specific  subdivision  schemes  was  presented  in  the  last 
two  chapters.  In  Chapter  3,  Catmull-Clark  subdivision  scheme,  a  representative  of 
the  approximating  subdivision  schemes,  is  embedded  in  a  physics-based  modeling  en- 
vironment. In  Chapter  4,  a  dynamic  framework  is  provided  for  the  modified  butterfly 
subdivision  scheme,  a  popular  interpolatory  subdivision  technique.  In  this  chapter, 
a  unified  approach  is  presented  for  embedding  any  subdivision  scheme  in  a  dynamic 
framework. 

5.1     Overview  of  the  Unified  Approach 

The  key  concept  of  the  unified  approach  is  the  ability  to  express  the  smooth 
limit  surface  generated  by  a  subdivision  scheme  as  a  collection  of  a  single  type  of 
finite  elements.  The  type  of  this  finite  element  depends  on  the  subdivision  scheme 
involved.  For  example,  the  limit  surface  generated  by  the  Catmull-Clark  subdivision 
scheme  can  be  expressed  as  a  collection  a  quadrilateral  finite  element  patches.  This 
approach  differs  form  the  one  taken  in  Chapter  3  where  the  limit  surface  was  a 
collection  of  normal  elements  (corresponding  to  the  quadrilateral  bicubic  B-spline 
patches  away  from  the  extraordinary  points)  and  special  elements  (corresponding  to 


83 


84 


the  n-sided  patches  near  extraordinary  vertices  of  degree  n).   This  concept  will  be 
further  elaborated  in  the  next  few  sections. 

The  unified  approach  adopts  two  different  strategies  depending  on  the  subdivi- 
sion scheme  involved.  One  strategy  is  for  the  approximating  schemes,  whose  limit 
surface  is  some  type  of  B-spline  surface  away  from  the  extraordinary  points,  and 
the  other  is  for  the  interpolatory  subdivision  schemes,  where  the  limit  surface  does 
not  have  any  analytical  expression.  It  may  be  recalled  that  recursive  subdivision  is 
a  method  for  smoothing  polyhedra,  and  hence  corresponding  to  each  n-sided  face 
in  the  control  mesh,  there  would  be  a  smooth  n-sided  patch  in  the  limit  surface 
in  general.  For  example,  the  control  mesh  (after  at  most  one  subdivision  step)  of 
the  Catmull-Clark  subdivision  scheme  consists  of  quadrilateral  faces,  which  when 
smoothed  through  subdivision  lead  to  quadrilateral  patches  in  the  limit  surface  (see 
Fig.5.1(b)).  In  the  case  of  approximating  subdivision  schemes,  the  faces  which  do 
not  contain  any  extraordinary  vertex  lead  to  spline  patches  with  analytic  expressions. 
On  the  other  hand,  faces  which  have  one  extraordinary  vertex,  lead  to  patches  whose 
analytic  expressions  are  difficult  to  determine.  Continuing  the  example  with  the 
Catmull-Clark  subdivision  scheme,  the  faces  without  extraordinary  vertices  lead  to 
quadrilateral  bicubic  B-spline  patches  in  the  limit,  whereas  the  faces  with  one  extraor- 
dinary vertex  lead  to  quadrilateral  patches  whose  analytic  expressions  are  difficult 
to  obtain.  In  Chapter  3,  n  such  adjacent  patches  corresponding  to  an  extraordinary 
vertex  of  degree  n  are  effectively  combined  together  to  obtain  an  expression  for  the  re- 
sulting n-sided  patch  (see  Fig.5.1(a)).  Recently,  Stam  [90,  91]  has  introduced  schemes 
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for  exactly  evaluating  the  limit  surface,  even  near  the  extraordinary  points,  for  ap- 
proximating subdivision  schemes.  Analytic  expressions  for  the  patches  resulting  from 
smoothing  faces  with  an  extraordinary  vertex  can  be  obtained  using  these  schemes. 
Therefore,  the  limit  surface  can  be  expressed  as  a  collection  of  a  single  patch  type, 
and  consequently  a  single  type  of  finite  element  can  be  used  to  provide  the  dynamic 
framework  for  the  approximating  subdivision  schemes.  This  unified  framework  for 
approximating  subdivision  schemes  is  fully  worked  out  in  this  chapter  for  Catmull- 
Clark  and  Loop  subdivision  schemes  and  a  general  outline  is  also  presented  on  how 
to  carry  out  the  same  strategy  for  other  approximating  subdivision  schemes. 

The  limit  surface  resulting  from  an  interpolators  subdivision  scheme  does  not 
have  an  analytic  expression  in  general,  and  hence  a  strategy  similar  to  the  one  adopted 
in  Chapter  4  needs  to  be  used.  An  outline  of  the  strategy  to  be  used  for  interpolatory 
subdivision  schemes  is  also  presented  in  this  chapter. 

5.2     Unified  Dynamic  Framework  for  Catmull-Clark  Subdivision  Surfaces 

A  systematic  formulation  along  with  the  implementation  details  of  the  unified 
dynamic  framework  for  Catmull-Clark  subdivision  surfaces  is  presented  in  this  sec- 
tion. The  key  difference  between  the  framework  developed  in  Chapter  3  and  the  one 
presented  here  is  the  representation  of  the  limit  surface,  which  leads  to  different  types 
of  finite  elements  used  for  developing  the  dynamic  framework.  This  is  illustrated  with 
a  schematic  diagram  in  Fig.5.1. 
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(a) 


(b) 


Figure  5.1.  A  control  mesh  with  an  extraordinary  vertex  of  degree  5  and  the  corre- 
sponding limit  surface  :  (a)  using  the  concepts  developed  in  Chapter  3,  where  the 
limit  surface  consists  of  quadrilateral  normal  elements  and  a  pentagonal  special  ele- 
ment; (b)  using  the  unified  concept  developed  in  this  chapter,  where  the  limit  surface 
consists  of  one  single  type  of  quadrilateral  finite  element. 
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Figure  5.2.  In  Catmull-Clark  subdivision,  each  non-boundary  quadrilateral  face  in 
the  control  mesh  has  a  corresponding  quadrilateral  patch  in  the  limit  surface  :  (a) 
control  mesh,  (b)  limit  surface. 


Following  the  concepts  developed  in  Chapter  3,  the  limit  surface  of  the  control 
mesh  shown  in  Fig.5.1,  consists  of  quadrilateral  bicubic  B-spline  patches  correspond- 
ing to  the  faces  marked  'n'  (faces  with  no  extraordinary  points),  and  a  pentagonal 
patch  corresponding  to  the  faces  marked  V  (faces  having  one  extraordinary  vertex 
of  degree  5)  (Fig.5.1(a)).  However,  in  this  section,  it  will  be  shown  that  the  entire 
limit  surface  can  be  expressed  as  a  collection  of  quadrilateral  patches  as  illustrated  in 
Fig.5.1(b)  using  the  schemes  proposed  by  Stam  [90].  Next,  the  local  parameterization 
of  the  limit  surface  is  described  in  order  to  develop  the  unified  dynamic  framework. 


5.2,1     Local  Parameterization 

As  mentioned  earlier,  the  control  mesh  (after  at  most  one  subdivision  step)  for 
the  Catmull-Clark  subdivision  scheme  consists  of  quadrilateral  faces  which  lead  to 
quadrilateral  patches  in  the  limit  surface.  For  the  sake  of  simplicity,  it  is  assumed 
that  each  face  has  at  most  one  extraordinary  vertex.  If  this  assumption  is  violated, 
then  one  more  subdivision  step  needs  to  be  performed  on  the  current  control  mesh 
in  order  to  obtain  a  new  control  mesh  on  which  the  following  analysis  can  be  carried 
out.  The  number  of  quadrilateral  patches  in  the  limit  surface  is  equal  to  the  number 
of  non-boundary  quadrilateral  faces  in  the  control  mesh  (Fig.5.2).  Therefore,  the 
smooth  limit  surface  s  can  be  expressed  as 


£s,,  (5-1) 

(=1 


where  n  is  the  number  of  non-boundary  faces  in  the  control  mesh  and  S;  is  the  smooth 
quadrilateral  patch  corresponding  to  the  /-th  non-boundary  quadrilateral  face  in  the 
control  mesh.  Each  of  these  quadrilateral  patches  can  be  parameterized  over  the 
corresponding  non-boundary  quadrilateral  face  in  the  control  mesh.  However,  since  a 
quadrilateral  face  can  easily  be  parameterized  over  a  [0, 1]  domain,  each  quadrilateral 
patch  is  locally  parameterized  over  the  domain  [0, 1]  . 

The  non-boundary  quadrilateral  faces  are  of  two  types  :  (a)  faces  having  no 
extraordinary  vertices  (dubbed  as  "regular"  faces  in  Section  3.2.1,  marked  as  k  in 
Fig.5.2(a))  and  (b)  faces  with  one  extraordinary  vertex  (dubbed  as  "irregular"  faces 
in  Section  3.2.2,  marked  as  Q  in  Fig.5.2(a)).  If  there  are  m  regular  and  n-m  irregular 
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faces,  then  Eqn.5.1  can  be  rewritten  as 


s  =  £s,+  £Sj,  (5.2) 


where  s;  is  the  quadrilateral  patch  corresponding  to  the  i-th  regular  face  and  s,  is 
the  quadrilateral  patch  corresponding  to  the  j'-th  irregular  face. 

The  quadrilateral  patch  in  the  limit  surface  corresponding  to  each  regular  face 
is  a  bicubic  B-spline  patch,  which  is  defined  over  [0, 1]  .  The  set  of  control  vertices 
defining  this  bicubic  B-spline  patch  can  be  obtained  using  the  adjacent  face  informa- 
tion (refer  to  Section  3.2.1  and  Fig.3.1  in  Chapter  3).  Therefore,  the  quadrilateral 
patches  in  the  smooth  limit  surface  corresponding  to  the  regular  faces  in  the  con- 
trol mesh  can  be  easily  expressed  analytically,  which  are  essentially  bicubic  B-spline 
patches  defined  by  16  control  vertices  over  a  [0,  l]2  domain.  The  analytic  expression 
for  the  quadrilateral  patch  corresponding  to  the  regular  face  i  is  given  by 

Si     =     J(,(U,  l))pi 

=    (Jt(u,u)Ai)p 

=    Jj(u,  v)p,  (5.3) 

where  0  <  u,v  <  1,  3t,(u,v)  is  the  collection  of  the  bicubic  B-spline  basis  functions, 
p,  is  the  concatenation  of  the  16  control  vertex  positions  defining  the  bicubic  B-spline 
patch,  Ai  is  a  selection  matrix  which  when  multiplied  with  p,  the  concatenation  of 
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all  the  control  vertex  positions  defining  the  smooth  limit  surface,  selects  the  corre- 
sponding set  of  control  vertices,  and  3i(u,v)  =  Jt(u,  v)A{. 

The  analytic  expression  of  the  quadrilateral  patches  corresponding  to  the  ir- 
regular faces  in  the  control  mesh  was  difficult  to  derive,  and  hence  an  alternative 
approach  was  taken  in  Mandal  et  al.  [57,  60]  and  Qin  et  al.  [76]  as  well  as  in  Chapter 
3.  However,  very  recently  an  efficient  scheme  for  evaluating  Catmull-Clark  subdivi- 
sion surfaces  at  arbitrary  parameter  values  was  proposed  by  Stam  [90],  The  proposed 
approach,  involving  eigen-analysis  of  the  subdivision  matrix,  leads  to  an  analytic  ex- 
pression of  a  quadrilateral  patch  parameterized  over  an  irregular  face.  Following  the 
scheme  developed  by  Stam  [90],  the  quadrilateral  patch  corresponding  to  the  irregular 
face  j  is  given  by 

Sj     =     J*  (".  1>)pj 

=    Jj(u,«)p,  (5.4) 

where  0  <  u, v  <  1  as  before.  J^t(u,v)  is  the  collection  of  basis  functions  for  the 
corresponding  quadrilateral  patch  in  the  smooth  limit  surface.  The  subscript  <4  is 
used  to  denote  the  fact  that  the  irregular  face  has  an  extraordinary  vertex  of  degree 
k.  The  detailed  derivation  and  the  analytic  expressions  of  these  basis  functions 
involving  the  eigenvalues  and  eigenvectors  of  the  subdivision  matrix  can  be  found  in 
Stam  [90].  The  rest  of  the  notations  used  in  Eqn.5.4  have  the  usual  meaning  :  pj  is 
the  concatenation  of  the  2fc  +  8  control  vertices  defining  the  quadrilateral  patch  in  the 
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Figure  5.3.  (a)  The  marked  16  control  vertices  define  the  shaded  quadrilateral  patch 
associated  with  the  shaded  regular  face  in  the  control  mesh;  (b)  the  marked  14  control 
vertices  define  the  shaded  quadrilateral  patch  associated  with  the  shaded  irregular 
face  in  the  control  mesh. 


limit  surface,  p  is  the  concatenation  of  all  the  control  vertex  positions  defining  the 
smooth  limit  surface,  Aj  is  a  selection  matrix  which  when  multiplied  with  p  selects 
the  corresponding  set  of  control  vertices,  and  3j(u,v)  =  3ik{u,v)Aj. 

It  may  be  noted  that  the  number  of  control  vertices  in  the  initial  mesh  defining 
a  quadrilateral  patch  in  the  smooth  limit  surface  is  Ik  +  8,  where  k  =  4  in  case 
the  associated  quadrilateral  face  in  the  control  mesh  is  regular,  or  k  =  degree  of  the 
extraordinary  vertex  if  the  associated  quadrilateral  face  is  irregular.  For  example,  the 
shaded  quadrilateral  patch  is  associated  with  the  shaded  regular  face  in  Fig.5.3(a), 
and  the  16  control  vertices  defining  this  patch  (which  is  actually  a  bicubic  B-spline 
patch)  are  marked.  Similarly,  the  shaded  quadrilateral  patch  is  associated  with  the 
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shaded  irregular  face  in  Fig. 5.3(b),  and  the  14  control  vertices  defining  this  patch  are 
highlighted.  Now  an  expression  of  the  smooth  limit  surface  can  be  easily  formulated. 
Using  Eqn.5.2,  5.3  and  5.4,  it  can  be  shown  that 


s 

=  EJ.p+EJiP 

t-i 

m 

=  (£J. 

+  E  J> 

=  Jp, 

(5.5) 


where  J  =  (££i  Ji  +  £?=™  Jj)-  Note  that  even  though  the  initial  mesh  serves  as  the 
parametric  domain  of  the  smooth  limit  surface,  each  quadrilateral  face  in  the  initial 
mesh  and  consequently  the  smooth  limit  surface  can  be  defined  over  a  [0, 1]    domain. 

5.2.2  Dynamics 

The  expression  of  the  limit  surface  as  given  by  Eqn.5.5  is  same  as  the  one  derived 
in  Chapter  3  (Eqn.3.15).  However,  the  contents  of  the  matrix  J  is  totally  different  in 
two  cases.  Nevertheless,  the  dynamics  of  the  limit  surface  can  be  formulated  using 
the  same  symbolic  expressions  (which,  when  actually  evaluated  will  lead  to  different 
expressions  depending  on  the  contents  of  the  matrix  J).  Therefore,  the  symbolic 
expression  of  the  motion  equation  remains  the  same  as  in  the  earlier  formulation  of 
the  dynamic  Catmull-Clark  subdivision  surface  model. 

5.2.3  Finite  Element  Implementation 

In  the  unified  dynamic  framework,  the  smooth  limit  surface  obtained  via  Catmull- 
Clark  subdivision  process  is  expressed  as  a  collection  of  quadrilateral  patches.  Each 
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quadrilateral  patch  in  the  limit  surface  is  considered  as  a  finite  element.  Therefore, 
the  limit  surface  under  the  unified  framework  consists  of  a  single  type  of  finite  el- 
ements instead  of  two  different  types  of  elements  as  in  Mandal  et  al.  [57,  60]  and 
Qin  et  al.  [76].  The  mass,  damping  and  stiffness  matrices  for  these  elements  along 
with  the  generalized  force  vector  can  be  evaluated  analytically,  provided  the  mate- 
rial properties  (like  mass,  damping,  rigidity  and  bending  distributions)  have  analytic 
expressions.  These  distribution  functions  are  assumed  to  be  constant  in  most  cases. 
Now,  the  expressions  of  the  mass,  damping  and  stiffness  matrices  for  a  quadrilateral 
element  which  is  a  bicubic  B-spline  can  be  written  as 

Me  =        /    niT3bdudv,  (5.6) 

Jo  Jo 

De  =  f  C  lJTbibdudv,  (5.7) 

Jo   Jo 

and 

Ke    =     f  /1(ou{(Jl).}r{(Jt).}  +  ffl22{(J»).}7'{(Jt).}  +  /8ii{(Jt)„}T{(Jt)M} 
Jo  Jo 

+/?i2{(J6Wr{(J(,)„»}  +  dMihUfUJuUVdudv  (5.8) 

respectively,  where  J&  is  the  bicubic  B-spline  basis  matrix,  /j(u,  v)  is  the  mass  density, 
•y(u,  v)  is  the  damping  density,  aa(u,  v)  and  f}ij(u,v)  are  the  tension  and  rigidity 
functions  respectively.  The  subscript  u  and  v  denote  partial  derivatives  with  respect 
to  u  and  v  respectively.  The  subscript  e  is  used  to  indicate  elemental  matrices  which 
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are  of  size  (16, 16).  The  expression  for  the  corresponding  generalized  force  vector  is 

fe=  f  [  3Tf(u,v,t)dudv.  (5.9) 

Jo  Jo 

The  mass,  damping  and  stiffness  matrices  for  the  quadrilateral  elements  which 
are  not  bicubic  B-splines  (corresponding  to  the  irregular  faces)  can  be  expressed 
analytically  as  well  simply  by  replacing  the  matrix  Jb  in  Eqn.5.6,  5.7  and  5.8  by 
the  matrix  3ik  (refer  to  Eqn.5.4),  where  k  denotes  the  degree  of  the  extraordinary 
vertex  associated  with  the  corresponding  irregular  face.  These  elemental  matrices 
are  of  size  (2k  +  8, 2k  +  8).  The  generalized  force  vector  for  these  elements  can  also 
be  determined  in  a  similar  fashion.  It  may  be  noted  that  the  limits  of  integration 
need  to  be  chosen  carefully  for  elemental  stiffness  matrices  as  the  second  derivative 
diverges  near  the  extraordinary  points  for  Catmull-Clark  subdivision  surfaces.  The 
elemental  matrices  combined  together,  either  explicitly  or  implicitly,  give  the  global 
matrices.  Similarly,  elemental  generalized  force  vectors  can  be  combined  together  to 
obtain  the  global  generalized  force  vector.  The  motion  equation  is  discretized  and 
solved  following  similar  techniques  described  in  Chapter  3. 

It  may  be  noted  that  even  though  an  analytical  expression  for  a  quadrilateral 
element  in  the  limit  surface  which  is  not  a  bicubic  B-spline  exists,  it  is  cumbersome 
to  actually  evaluate  the  elemental  matrix  expressions.  Numerical  integration  using 
Gaussian  quadrature  may  be  used  to  obtain  approximations  of  these  elemental  ma- 
trices.  However,  in  this  dissertation,  an  approach  similar  to  the  dynamic  modified 
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butterfly  subdivision  model  presented  in  Chapter  4  is  taken  for  implementation  pur- 
poses. An  approximation  of  the  smooth  limit  surface  is  obtained  by  subdividing 
the  initial  control  mesh  j  times,  and  a  spring-mass  system  is  developed  on  this  j-th 
approximation  level  in  a  similar  fashion  as  in  Section  4.3.  The  physical  matrices 
of  this  system  are  used  as  approximations  to  the  actual  physical  matrices.  These 
approximations  have  been  found  to  be  very  efficient  for  implementation  purposes. 

5.3     Unified  Dynamic  Framework  for  Loop  Subdivision  Surfaces 

Loop's  subdivision  scheme  starts  with  a  triangular  control  mesh  and  generates  a 
smooth  surface  with  triangular  patches  in  the  limit.  It  is  an  approximating  subdivi- 
sion scheme  which  generalizes  recursive  subdivision  schemes  for  obtaining  C2  quartic 
triangular  B-spline  patches  in  a  regular  setting.  In  each  step  of  Loop  subdivision,  each 
(non-boundary)  triangular  face  is  refined  into  4  triangular  faces  using  the  following 
rules  : 

•  For  each  (non-boundary)  vertex  V  of  degree  n,  a  new  vertex  point  is  introduced. 
The  position  of  this  newly  introduced  vertex  point  is  given  by  °  „Z]+n "+v" , 
where  v  is  the  position  vector  of  vertex  V,  Vj, . . . ,  vn  are  the  vertex  positions 
of  the  n  vertices  connected  to  vertex  V,  a(n)  =  "^ff1^  and  (3(n)  =  |  - 

(3+2cos(27r/n))2 
64 

•  For  each  (non-boundary)  edge  E,  a  new  edge  point  is  introduced.  Let  E  be  the 
connecting  edge  between  vertices  V\  and  V2,  and  is  shared  by  faces  F\  and  Fi. 
If  Fi  and  F2  have  vertices  VFl  and  Vp2  respectively  (apart  from  Vi  and  V'2),  then 
the  position  of  the  newly  introduced  edge  point  is  given  by  8VF'    "F2 , 
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(a) 


(b) 


Figure  5.4.   (a)  The  control  polygon  with  triangular  faces;  (b)  mesh  obtained  after 
one  subdivision  step  using  Loop's  subdivision  rules. 


where  Vi,  v2,  v^1  and  vp,  are  the  position  vector  of  the  vertex  Vi,  Vi,  Vf,  and 
Vf2  respectively. 

•  New  edges  are  formed  by  connecting  each  new  vertex  point  to  the  new  edge 
points  corresponding  to  the  edges  incident  on  the  old  vertex,  and  by  connecting 
each  new  edge  point  to  the  new  edge  points  of  the  other  edges  in  the  two  faces 
which  shared  the  original  edge. 

•  New  faces  are  defined  as  faces  enclosed  by  the  new  edges. 

Examples  of  refining  an  initial  mesh  using  Loop's  subdivision  rules  are  shown 
in  Fig.5.4  and  5.5.  These  subdivision  rules  ensure  tangent  plane  continuity  of  the 
limit  surface  even  in  a  irregular  setting,  i.e.,  when  the  triangular  control  mesh  has 
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extraordinary  vertices  whose  degree  is  not  equal  to  6.  The  limit  surface  shown  in 
Fig.2.1(d)  is  actually  obtained  by  subdividing  the  control  mesh  shown  in  Fig.2.1(a) 
using  Loop  subdivision  rules.  A  detailed  discussion  on  how  to  obtain  positions  and 
normals  in  the  smooth  limit  surface  generated  by  the  Loop  subdivision  scheme  can 
be  found  in  Section  4.2  of  Hoppe's  thesis  [39]. 

5.3.1     Local  Parameterization 

The  limit  surface  obtained  via  Loop's  subdivision  scheme  can  be  locally  pa- 
rameterized easily.  This  local  parameterization  scheme  is  very  similar  in  nature  to 
the  one  described  for  Catmull-Clark  subdivision  scheme  in  Section  5.2.1.  For  Loop's 
scheme,  the  smooth  limit  surface  consists  of  triangular  patches  and  the  number  of 
these  triangular  patches  is  same  as  the  number  of  non-boundary  triangular  faces  in 
the  control  mesh.  Therefore,  each  of  the  triangular  patch  in  the  limit  can  be  locally 
parameterized  over  the  corresponding  triangular  face  in  the  control  mesh.  In  may 
be  noted  that  each  triangular  face  in  the  control  mesh  can  be  parameterized  over  a 
triangular  domain  whose  vertices  are  located  at  (0,0),  (0,1)  and  (1,0),  and  hence 
each  triangular  patch  and  consequently  the  smooth  limit  surface  can  be  defined  over 
this  domain  (refer  Fig.5.6). 

The  triangular  patches  in  the  smooth  limit  surface  are  of  two  types.  For  a  non- 
boundary  triangular  face  in  the  control  mesh  with  no  extraordinary  vertices  (i.e.,  with 
three  vertices  of  degree  6),  the  corresponding  triangular  patch  in  the  limit  surface  is 
a  particular  type  of  triangular  B-spline  (the  three-direction  quartic  box  spline)  whose 
analytic  expression  is  easy  to  obtain.  This  triangular  B-spline  patch  is  controlled  by 
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control  mesh 


limit  surface 


Figure  5.5.  An  initial  mesh  and  the  corresponding  limit  surface  obtained  using  Loop's 
subdivision  rules.  The  domains  of  the  shaded  triangular  patches  in  the  limit  surface 
are  the  corresponding  triangular  faces  in  the  initial  mesh.  The  encircled  vertices  are 
the  control  vertices  for  the  corresponding  triangular  patch  in  the  limit  surface. 
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(0,1 


(0,0)         (1,0) 


initial  rnesh 


4  \ 


limit  surface 


Figure  5.6.  Each  triangular  patch  in  the  limit  surface  can  be  associated  with  a  non- 
boundary  triangular  face  in  the  initial  mesh,  which  in  turn  can  be  parameterized  over 
a  triangle  with  vertices  at  (0,0),  (1,0)  and  (0, 1). 
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12  vertices  as  shown  in  Fig.5.5  (the  set  of  enclosed  vertices  in  the  left  hand  side).  The 
triangular  patch  in  the  limit  surface  corresponding  to  a  non-boundary  triangular  face 
in  the  control  mesh  with  one  extraordinary  vertex  can  also  be  expressed  analytically 
using  the  schemes  proposed  by  Stam  [91].  This  triangular  patch  is  controlled  by  n  +  6 
vertices  in  the  control  mesh  where  n  is  the  degree  of  the  extraordinary  vertex.  The 
set  of  control  vertices  for  a  triangular  patch  of  the  later  type  is  shown  in  the  right 
hand  side  of  Fig.5.5.  Therefore,  each  triangular  patch  in  the  limit  surface  can  be 
expressed  analytically,  and  an  expression  of  the  limit  surface  similar  to  Eqn.5.5  can 
be  obtained. 

5.3.2  Dynamics 

Once  an  expression  for  the  limit  surface  obtained  via  Loop's  subdivision  is  ob- 
tained, the  dynamic  model  can  be  developed  following  an  exactly  similar  procedure 
described  for  Catmull-Clark  subdivision  scheme  in  Section  5.2.2.  The  motion  equa- 
tion of  the  dynamic  Loop  subdivision  model  can  also  be  derived  in  a  similar  fashion. 

5.3.3  Finite  Element  Implementation 

The  implementation  of  the  dynamic  framework  for  Loop  subdivision  scheme  us- 
ing the  unified  approach  treats  each  triangular  patch  in  the  limit  surface  as  a  finite 
element.  Each  triangular  patch  has  an  analytic  expression,  and  hence  the  elemental 
physical  matrices  and  the  generalized  force  vector  can  be  derived  analytically.  The 
derivation  of  an  exact  expression  for  elemental  matrices  is  cumbersome  for  the  tri- 
angular patches  corresponding  to  the  triangular  faces  with  an  extraordinary  vertex, 
and  numerical  integration  using  Gaussian  quadrature  may  be  used  for  deriving  an 
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approximation.  However,  a  practical  alternative  for  implementation  is  to  subdivide 
the  control  mesh  j  times  using  Loop's  subdivision  rules,  and  to  build  a  spring-mass 
system  on  this  j-th  level  approximation  as  has  been  done  for  the  dynamic  modified 
butterfly  subdivision  model  in  Section  4.3.  The  physical  matrices  of  this  spring-mass 
system  provide  an  approximation  of  the  original  physical  matrices,  and  it  works  well 
in  practice. 

5.4     A  General  Outline  of  the  Framework  for  Approximating  Subdivision  Schemes 

The  unified  approach  presented  to  provide  a  dynamic  framework  for  Catmull- 
Clark  and  Loop  subdivision  schemes  can  be  generalized  easily.  This  approach  involves 
three  steps  which  are  as  follows  : 

1.  The  limit  surface  obtained  via  an  approximating  subdivision  scheme  can  be 
expressed  as  a  collection  of  smooth  patches  which  can  be  locally  parameterized 
over  a  corresponding  face  in  the  control  mesh.  Each  patch  is  n-sided  if  it 
is  locally  parameterized  over  a  n-sided  face.  Analytic  expressions  for  each  of 
these  patches  can  be  derived  even  in  the  presence  of  extraordinary  vertices  in 
the  control  mesh,  and  hence  an  expression  of  the  limit  surface  can  be  obtained. 

2.  Once  an  expression  of  the  limit  surface  is  obtained,  the  dynamic  framework  can 
be  developed  by  treating  control  vertex  positions  as  a  function  of  time.  The 
corresponding  motion  equation  can  be  derived. 
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3.  Each  patch  in  the  limit  surface  is  treated  as  a  finite  element  in  implementation. 
The  elemental  mass,  damping  and  stiffness  matrices  along  with  the  general- 
ized force  vector  can  be  obtained  by  either  analytic  or  numerical  integration. 
Alternatively,  the  control  mesh  can  be  subdivided  j  times  to  obtain  an  approxi- 
mation of  the  smooth  limit  surface,  and  a  spring-mass  system  can  be  developed 
on  this  approximation  mesh.  The  physical  matrices  of  this  system  provide  an 
approximation  to  the  original  physical  matrices  and  works  well  in  practice. 

5.5     A  General  Outline  of  the  Framework  for  Interpolatory  Subdivision  Schemes 

Most  of  the  interpolatory  subdivision  schemes  are  obtained  by  modifying  the 
butterfly  subdivision  scheme  [30].  Therefore,  the  framework  developed  for  the  modi- 
fied butterfly  subdivision  scheme  in  Chapter  4  is  pretty  general  and  can  be  used  for 
other  interpolatory  subdivision  schemes  as  well.  The  only  difference  is  that  the  basis 
functions  as  well  as  the  set  of  control  vertices  for  a  patch  in  the  limit  surface  depend 
on  the  chosen  interpolatory  subdivision  rules.  It  may  also  be  noted  that  unlike  the 
approximating  schemes,  the  physical  matrices  can  not  be  obtained  analytically  as 
the  basis  functions  corresponding  to  interpolatory  subdivision  schemes  do  not  have 
any  analytic  expressions  in  general.  Even  though  these  matrices  can  be  obtained  via 
numerical  integration,  the  spring-mass  system  developed  in  Chapter  4  is  preferred 
for  implementation  purposes  due  to  efficiency  reasons. 


CHAPTER  6 
MULTIRESOLUTION  DYNAMICS 


The  dynamic  framework  developed  for  subdivision  surfaces  allows  the  users  to 
directly  manipulate  the  smooth  limit  surface  by  applying  synthesized  forces.  In  any 
subdivision  scheme,  an  initial  (control)  mesh  is  refined  using  specific  subdivision 
rules  to  obtain  a  smooth  surface  in  the  limit.  To  develop  the  dynamic  framework, 
this  smooth  limit  surface  is  parameterized  over  the  domain  defined  by  the  initial  mesh 
and  hence,  the  limit  surface  can  be  expressed  as  a  function  of  the  initial  mesh.  The 
users  apply  synthesized  forces  on  the  smooth  limit  surface  which  is  transferred  back 
on  the  initial  mesh  via  a  transformation  matrix.  The  initial  mesh  evolves  over  time  in 
response  to  the  applied  forces,  and  consequently  the  smooth  limit  surface  deforms  to 
obtain  an  equilibrium  position  where  all  forces  are  balanced.  As  noted  earlier,  even 
though  various  types  of  forces  can  be  applied  on  the  limit  surface  to  obtain  a  desired 
effect,  the  possible  shapes  that  can  be  obtained  via  evolution  is  directly  related  to 
the  number  of  control  vertices  present  in  the  initial  mesh. 

The  concept  of  multilevel  dynamics  was  introduced  earlier  in  order  to  have  vary- 
ing number  of  control  vertices  in  the  initial  mesh,  and  a  wide  range  of  shapes  with 
arbitrary  topology  can  be  modeled  within  the  dynamic  framework.    The  reader  is 
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Figure  6.1.  Schematic  block  diagram  of  the  multilevel  dynamics  approach. 
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referred  to  Section  3.2.5  and  Section  3.3.6  of  this  dissertation  for  the  details  of  multi- 
level dynamics  in  the  context  of  dynamic  Catmull-Clark  surfaces,  and  to  Section  4.2.3 
for  the  dynamic  modified  butterfly  subdivision  scheme.  The  concept  of  multilevel  dy- 
namics is  illustrated  with  block  diagrams  in  Fig.6.1.  This  multilevel  dynamics  can 
be  considered  as  a  subdivision  surface-based  coarse-to-fine  modeling  framework  and 
can  be  summarized  as  follows.  The  user  starts  with  a  smooth  limit  surface  which  has 
a  very  simple  initial  mesh  with  few  control  vertices.  This  smooth  surface  is  deformed 
by  applying  synthesized  forces  so  that  the  limit  surface  at  equilibrium  would  look 
like  the  final  desired  shape  but  without  the  details  in  it.  Then,  the  number  of  control 
vertices  in  the  initial  mesh  is  increased  by  replacing  the  initial  mesh  at  equilibrium 
by  another  one  obtained  via  one  subdivision  step  on  the  old  initial  mesh,  and  the 
old  initial  mesh  is  discarded.  Both  of  these  old  and  new  initial  mesh  represent  the 
same  limit  surface,  but  the  latter  has  more  control  vertices.  Now,  this  new  initial 
mesh,  and  consequently  the  smooth  limit  surface,  can  be  deformed  to  obtain  a  shape 
at  equilibrium  which  has  more  details.  This  process  is  repeated  till  the  deformed 
smooth  surface  has  all  the  desired  details. 

The  limitation  of  the  multilevel  dynamics  approach  is  that  once  a  finer  resolution 
is  chosen,  the  modeler  can  not  opt  for  a  lower  resolution  as  the  dynamics  is  defined  on 
the  current  initial  mesh.  There  is  no  way  of  having  the  dynamics  defined  on  a  simpler 
initial  mesh,  at  the  same  time  preserving  the  already  obtained  details  in  the  model 
within  the  multilevel  dynamics  approach.  This  is  potentially  a  severe  restriction.  For 
example,  let  us  imagine  the  scenario  depicted  in  Fig.6.2.  The  modeler  has  recovered 
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Figure  6.2.  Multilevel  dynamics  vs.  multiresolution  dynamics. 
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a  long  straight  pipe-like  shape  with  lots  of  detail  by  synthesized  force  application 
on  a  smooth  limit  surface.  This  smooth  limit  surface  initially  had  a  simple  initial 
mesh  with  few  vertices,  but  the  current  initial  mesh  after  many  model  subdivision 
steps  has  large  number  of  control  vertices.  Now,  if  the  modeler  wants  to  bend  this 
straight  pipe-like  shape,  the  synthesized  forces  applied  on  the  smooth  limit  surface 
need  to  deform  an  initial  mesh  with  a  large  number  of  control  vertices.  It  would  have 
been  much  easier  and  faster  if  a  lower  resolution  version  of  the  current  initial  mesh 
could  be  deformed  with  the  details  being  added  back  when  the  shape  has  bent  to 
the  desired  extent  in  response  to  the  applied  synthesized  forces.  Also,  it  has  been 
pointed  out  in  Section  2.5  that  the  dynamic  framework  can  integrate  shape  recovery 
and  shape  modeling  in  a  unified  framework  where  the  recovered  shape  can  easily  be 
edited  as  well.  However,  if  the  modeler  wants  a  global  change  in  the  recovered  shape, 
a  lower  resolution  control  mesh  (than  the  one  obtained  via  shape  recovery  process) 
may  be  required,  and  the  details  need  to  be  preserved  at  the  same  time. 

A  multiresolution  representation  of  the  initial  mesh  is  proposed  in  this  chapter 
to  solve  the  above-mentioned  problem.  At  equilibrium,  an  initial  mesh  is  subdivided 
to  obtain  an  initial  mesh  with  more  control  vertices  (degrees  of  freedom).  Instead 
of  developing  the  dynamics  on  this  newly  obtained  initial  mesh  as  in  the  multilevel 
dynamics  approach,  the  dynamics  is  developed  on  a  representation  which  views  the 
newly  obtained  initial  mesh  as  the  previous  initial  mesh  plus  detail  (wavelet)  co- 
efficients needed  to  get  the  current  initial  mesh  (see  Fig. 6. 3).  This  is  essentially 
developing  the  dynamics  on  a  multiresolution  representation  of  the  evolving  initial 
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degrees  of  freedom  =  3 
( 3  control  vertices) 


degrees  of  freedom  =  6 
( 6  control  vertices) 


(a)  Multilevel  approach 


degrees  of  freedom  =  3 
( 3  control  vertices) 


degrees  of  freedom  =  6 

(  3  control  vertices,  3  detail  coefficients) 


(b)  Multiresolution  approach 


Figure  6.3.    Representation  of  the  degrees  of  freedom  in  multilevel  dynamics  and 
multiresolution  dynamics  approach. 
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mesh  which  has  undergone  model  subdivision  step(s)  to  increase  the  degrees  of  free- 
dom. It  may  be  noted  that  both  the  multilevel  and  the  multiresolution  approach 
have  same  degrees  of  freedom  to  represent  the  initial  mesh  obtained  after  a  model 
subdivision  step,  but  the  degrees  of  freedom  for  the  multilevel  dynamics  approach  are 
the  control  vertex  positions  in  the  newly  obtained  mesh  whereas,  the  degrees  of  free- 
dom for  the  current  approach  are  the  control  vertex  positions  in  the  previous  initial 
mesh  and  the  detailed  coefficients  necessary  to  obtain  the  current  initial  mesh  from 
the  previous  one.  The  multiresolution  approach  has  the  benefit  that  the  modeler  can 
opt  for  a  lower  resolution  initial  mesh  later,  on  which  the  dynamics  can  be  defined, 
at  the  same  time  preserving  the  detailed  coefficients  of  the  initial  mesh  on  which  the 
dynamics  is  currently  defined.  This  concept  will  be  detailed  further  in  later  sections, 
but  first  an  overview  of  the  multiresolution  analysis  is  provided  so  that  the  reader 
can  easily  follow  the  rest  of  the  material  presented  in  this  chapter. 

6.1     Overview  of  Multiresolution  Analysis  and  Wavelets 

The  basic  idea  of  multiresolution  analysis  is  to  decompose  a  complicated  function 
into  a  lower  resolution  version  along  with  a  set  of  detailed  coefficients  (also  known 
as  wavelet  coefficients)  necessary  to  recover  the  original  function.  Multiresolution 
analysis  and  wavelets  have  widespread  applications  in  diverse  areas  like  signal  analysis 
[55,  81],  image  processing  [26,  56],  physics  [2],  numerical  analysis  [7],  computer  vision 
[48]  and  computer  graphics  [6,  13,  16,  17,  31,  32,  35,  45,  50,  53,  67,  68,  70,  85, 
86,  104,  110].  The  specific  areas  in  computer  graphics  that  involve  multiresolution 
analysis  and  wavelets  include,  but  not  limited  to,  global  illumination  [16,  17,  35,  85], 
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curve,  mesh  and  image  editing  [6,  32,  110],  video  processing  [31],  surface  viewing  [13], 
surface  interpolation  [70],  creating  surfaces  from  contours  [67,  68],  animation  control 
[50],  polyhedral  compression  [53],  modeling  [45,  104]  and  representing  functions  on 
a  sphere  [86].  Various  applications  of  wavelets  in  computer  graphics  are  detailed  in 
Stollnitz  et  al.  [92].  The  idea  of  multiresolution  analysis  is  explained  next  with  a 
simple  example. 

Let  p>  be  a  vector  containing  functional  values  of  some  arbitrary  function  at  n 
discrete  points.  A  lower  resolution  representation  of  the  function  using  m  (m  <  n) 
functional  values  can  be  obtained  from  p1  by  doing  certain  weighted  averages  of  the 
n  functional  values  stored  in  p>.  This  is  essentially  a  low  pass  filtering  followed  by 
down-sampling,  and  the  entire  operation  can  be  expressed  as 

p>->  =  AV,  (6-1) 

where  p7-1  is  the  m-valued  vector  representing  a  lower  resolution  version  and  the 
matrix  A3  is  of  size  (m,  n)  and  implements  the  low  pass  filter  along  with  the  down- 
sampler. 

The  lower  resolution  p*-1  has  fewer  functional  values  than  p>,  and  hence  some 
amount  of  detail  present  in  p1'  is  lost  due  to  the  low  pass  filtering  process.  If  the  low 
pass  filter  AJ  is  chosen  appropriately,  it  is  possible  to  capture  the  lost  detail  by  high 
pass  filtering  followed  by  a  down-sampling  of  the  original  functional  values  stored  in 
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p*.  This  operation  can  be  expressed  as 

tf-1  =  BV,  (6-2) 

where  q-*-1  is  a  (n  -  m)-valued  vector  representing  detailed  coefficients  necessary  to 
recover  p7  from  p*_1.  The  matrix  BJ  is  of  size  (n  -  m,n)  and  implements  the  high 
pass  filter  along  with  the  down-sampler.  The  high  pass  filter  BJ  is  related  to  the  low 
pass  filter  AJ.  These  two  filters,  AJ  and  BJ,  are  known  as  analysis  filters  and  the 
process  of  splitting  a  high  resolution  vector  p*  into  a  low  resolution  vector  pJ-1  and 
a  vector  qJ_1  storing  detail  coefficients  is  known  as  decomposition. 

The  original  vector  p*  can  be  reconstructed  from  the  vectors  p*1  and  qJ_1  if 
the  synthesis  filters  are  chosen  correctly.  This  is  called  the  reconstruction  process, 
which  can  be  expressed  as 

p>  =  TV-1  +  Qty~\  (6-3) 

where  Pj  is  a  refinement  matrix  of  size  (n,  m)  and  QJ  is  a  perturbation  matrix  of 
size  (n,m  -  n).  The  refinement  matrix  Pj  encodes  the  rules  for  obtaining  a  vector 
of  larger  size  from  a  given  vector  and  the  perturbation  matrix  QJ  encodes  how  to 
obtain  a  perturbation  of  (pJ  -  P^p*-1)  from  the  given  detailed  coefficient  vector  qJ_I. 
The  matrices  PJ  and  QJ  are  collectively  called  synthesis  filters,  and  are  related  to 
the  analysis  filters. 
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qJ  qJ 

Figure  6.4.  The  filter  bank. 


The  process  of  decomposition  can  be  continued  recursively.  For  example,  p*  is 
decomposed  into  a  lower  resolution  part  p>  '  and  a  detail  part  q»~  .  Then,  p*_1 
can  be  expressed  as  a  lower  resolution  part  p>-2  and  another  detail  part  qJ_2.  This 
recursive  process  is  depicted  in  Fig. 6. 4  and  is  known  as  filter  bank.  The  original  vector 
is  decomposed  into  a  hierarchy  of  lower  resolution  parts  p*-1 ,  p*-2, . . . ,  p°  and  detail 
parts  qJ_1,  qJ_2, . . . ,  q°.  The  original  vector  p*  can  be  recovered  from  the  sequence 
p°,q°,ql, .  • .  ,ij'_2,q,'_1.  This  sequence  has  the  same  size  as  that  of  the  original 
vector  and  is  known  as  the  wavelet  transform  of  the  original  vector. 

The  analysis  and  synthesis  filters  are  designed  in  such  a  fashion  that  the  low 
resolution  versions  are  good  approximations  of  the  original  function  in  a  least  square 
sense.  The  decomposition  and  reconstruction  processes  should  have  time  complexity 
0(n),  where  n  is  the  size  of  the  vector  being  decomposed  or  reconstructed.  A  de- 
tail (wavelet)  coefficient  value  should  also  be  related  to  the  error  introduced  in  the 
approximation  when  that  particular  coefficient  is  set  to  zero. 
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Multiresolution  analysis  is  a  framework  for  developing  these  analysis  and  syn- 
thesis filters.  It  involves  derivation  of  a  sequence  of  nested  linear  spaces  V'c  Vc 
V2  C  . . .  such  that  the  resolution  of  the  functions  contained  in  some  space  VJ  in- 
creases with  increasing  j.  Also,  there  exists  an  orthogonal  complement  space  W'~l 
for  each  V'~l,  j  =  1,  2, . . .  such  that  V'~l  ©  W'~l  =  Vj,  i.e.,  the  linear  space  VJ_1 
and  its  orthogonal  complement  space  W'~l  together  span  the  linear  space  VK  The 
inner  product  between  any  two  functions  at  some  level  j,  j  =  0, 1, 2, . . .  needs  to  be 
defined  in  order  to  derive  the  orthogonal  complement  space. 

Let  <jr{,  </4,  ■ . .,  <j^n  be  a  set  of  basis  functions  spanning  the  linear  space  VK 
These  basis  functions  are  also  called  scaling  functions  for  level  j.  Now,  VJ_1  c  VJ 
and  hence  the  basis  functions  $_1,  <^~  ,  . . .,  4>mXi  spanning  the  linear  space  V>~x, 
can  be  expressed  as  a  linear  combination  of  the  basis  functions  than  span  the  linear 
space  VK  This  refinability  of  scaling  functions  is  used  to  construct  the  refinement 
matrix  PJ,  which  should  satisfy  the  expression 

where  &1  =  [^j"1,^"1, . . . ,  (fa1],  *J  =  [<#,<&  . . .  ,(fa\,  and  P>  is  the  (n,m)  re- 
finement matrix  introduced  earlier  (Eqn.6.3).  Similarly,  let  J^i"1,  j/>2_1>  ■  ■  ■!  ^n-m  be 
a  set  of  basis  functions  that  span  the  linear  space  W^~l.  These  basis  functions,  also 
known  as  wavelets  at  level  j  —  1,  together  with  the  basis  functions  of  linear  space 
Vi~l  form  a  basis  for  the  linear  space  VK  Next,  the  perturbation  matrix  Q^  can  be 
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constructed  such  that  it  satisfies  the  expression 

&-1  =  &Qi,  (6.5) 

where  *J_1  =  [${  ,$j  ,  •  ■  • ,  i4i-m]t  ar,d  Q-*  is  the  (n,  n  -  m)  perturbation  matrix 
introduced  in  Eqn.6.3.  A  more  compact  expression  may  be  obtained  by  combining 
Eqn.6.4  and  6.5,  which  can  be  written  as 

[&-1  |  ^J"1]       =       &    \p*  |  Q>].  (6.6) 

The  analysis  filters  AJ  and  BJ  should  satisfy  the  inverse  relation 


Both  the  matrices  [Pj  |  QJ]  and    gj    are  of  size  (n,n)  and  satisfy  the  relation 

P>IQ>]     =      ff 


(6.7) 


(6.8) 


6.2     Multiresolution  Analysis  for  Surfaces  of  Arbitrary  Topology 

Classically,  the  multiresolution  analysis  was  developed  on  infinite  domains  such 
as  real  line  5R  and  plane  S{2  [22,  23].  Infinite  domains  are  spatially  invariant  and 
hence  a  single  basis  function  at  the  coarsest  resolution  can  be  dilated  and  translated 
to  form  the  basis  functions  at  higher  resolutions.  However,  many  applications  need 
multiresolution  analysis  on  bounded  intervals  and  techniques  for  imposing  boundary 


115 

constrains  have  also  been  derived  [19,  79].  The  idea  of  generalizing  multiresolution 
analysis  on  arbitrary  manifolds  was  first  introduced  by  Lounsbery  et  al.  [52,  53],  and 
was  further  improved  by  Stollnitz  et  al.  [92]  using  the  concepts  of  lifting  scheme  [93], 
In  this  dissertation,  the  multiresolution  schemes  derived  by  Stollnitz  et  al.  [92]  for 
surfaces  of  arbitrary  topology  will  be  adopted  to  develop  the  multiresolution  dynamic 
framework  for  subdivision  surfaces. 

Multiresolution  analysis  using  wavelets  can  be  categorized  into  three  types  - 
orthogonal  [22,  55],  semiorthogonal  [18]  and  biorthogonal  [20,  93].  A  multiresolution 
scheme  has  orthogonal  wavelets  if  the  wavelets  are  orthogonal  to  one  another,  and 
every  wavelet  is  orthogonal  to  every  coarser  scaling  function.  It  is  difficult  to  con- 
struct orthogonal  wavelets  that  have  local  support  (i.e.,  non-zero  over  a  small  section 
of  the  domain).  Locally  supported  wavelets  are  important  as  they  lead  to  sparse 
synthesis  and  analysis  filters  which  in  turn  lead  to  linear  time  complexity  in  vector 
size.  Wavelets  with  local  support  can  be  constructed  by  relaxing  the  orthogonality 
constraint.  If  the  wavelets  are  orthogonal  to  each  other  within  one  resolution  but 
not  across  different  resolutions,  then  they  are  called  semiorthogonal  wavelets.  In 
a  biorthogonal  setting,  the  orthogonality  constraint  is  dropped  where  W^~x  is  some 
(and  not  orthogonal)  complement  of  V'~l  in  VJ.  The  wavelets  on  arbitrary  manifolds 
developed  by  Stollnitz  et  al.  [92]  and  used  in  this  dissertation  are  biorthogonal.  The 
multiresolution  analysis  on  arbitrary  manifolds  involve  (1)  developing  nested  spaces 
using  subdivision  schemes,  (2)  selecting  an  inner  product  on  arbitrary  manifold  and 
(3)  constructing  biorthogonal  wavelets,  which  are  discussed  in  the  rest  of  this  section. 
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It  may  be  noted  that  there  is  a  slight  deviation  from  the  standard  notational 
convention  of  writing  vectors  in  bold  face  lower  case  and  matrices  in  bold  face  upper 
case  letters  in  a  particular  situation  in  the  rest  of  this  chapter.  The  collection  of 
control  point  positions  are  assumed  to  be  a  three-column  matrix  (one  for  each  of 
the  x,  y  and  z  components)  instead  of  a  vector  whose  size  is  thrice  the  number  of 
control  points.  Consequently,  the  three-row  matrix  representing  the  collection  of 
basis  functions  becomes  a  row  vector.  This  approach  is  more  natural  for  developing 
the  wavelet  expressions,  whereas  the  opposite  was  true  for  the  theory  developed  in 
earlier  chapters.  However,  bold  face  lower  case  letters  are  used  to  represent  the  three- 
column  matrix  of  control  vertex  positions  and  bold  face  upper  case  letters  are  used  to 
represent  the  vector  of  basis  functions,  so  that  the  reader  does  not  get  confused  with 
a  new  notation  for  representing  the  same  thing  in  a  different  format.  The  notational 
purity  is  given  away  to  a  small  extent  for  the  sake  of  clarity.  All  the  expressions  can  be 
reformulated  by  merging  the  columns  of  the  three-column  matrix  into  a  single  vector, 
but  the  resulting  expressions  will  unnecessarily  obscure  the  intuitiveness  contained 
in  the  current  set  of  expressions. 

6.2.1     Nested  Spaces  through  Subdivision 

In  previous  chapters,  it  has  been  shown  that  the  smooth  limit  surface  generated 
by  a  recursive  subdivision  scheme  in  general  can  be  written  as 


s(x)  =  J(x)p,  (6.9) 
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where  s  is  the  smooth  limit  surface,  J  is  a  collection  of  basis  functions  spanning 
the  space  defined  by  the  control  mesh  S,  p  is  the  vertex  positions  in  the  control 
mesh  and  x  €  S.  A  general  outline  on  how  to  obtain  this  expression  for  any  type  of 
subdivision  scheme  has  been  provided  in  Section  5.4  and  5.5  of  this  dissertation  and 
specific  examples  can  be  found  in  Eqn.3.15  and  Eqn.5.5  for  Catmull-Clark  subdivision 
scheme  and  in  Eqn.4.11  for  the  modified  butterfly  scheme. 

Let  the  control  mesh  defining  the  smooth  limit  surface  be  5°.  Renaming  the 
control  vertex  vector  p  to  p°  and  the  basis  (scaling)  function  matrix  J  to  J°,  Eqn.6.9 
can  be  rewritten  as 

s(x)  =  J°(x)p°,  (6.10) 

where  x  £  S°.  Let  S1  be  a  mesh  obtained  by  subdividing  5°  once  using  the  sub- 
division rules.  If  the  control  vertices  of  the  mesh  S1  are  collected  in  a  vector  p1, 
then 

p'  =  PV,  (6.11) 

where  P1  is  the  subdivision  (refinement)  matrix.  However,  the  control  mesh  Sl  has 
the  same  limit  surface  (as  it  is  obtained  via  one  subdivision  step  on  the  control  mesh 
S°)  and  hence  the  same  limit  surface  can  be  written  as 

s(x)    =    j'Wp1 

=    J'MPV,  (6.12) 
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where  J1  is  a  collection  of  basis  (scaling)  functions  spanning  the  space  defined  by  S1. 
From  Eqn.6.10  and  Eqn.6.12,  it  is  obvious  that 

J°(x)  =  J1(x)P1,  (6.13) 

which  essentially  means  that  the  scaling  functions  at  subdivision  level  0  can  be  ex- 
pressed as  a  linear  combination  of  the  scaling  functions  at  subdivision  level  1,  i.e., 
V°(S°)  C  V1(S°).  This  can  be  generalized  for  other  subdivision  levels  as  well,  and 
hence  subdivision  naturally  leads  to  a  nested  space  sequence  V°(S°)  C  V1(S°)  C 
V2(S°)  C  . . .  which  will  be  utilized  for  multiresolution  analysis  on  arbitrary  mani- 
folds. 

6.2.2     Inner  Product 

The  inner  product  proposed  in  Stollnitz  et  al.  [92]  is  used  in  this  dissertation. 
If  /„  and  /&  are  two  functions  defined  on  the  domain  5°,  then  their  inner  product 
<  fa,fb>  is  given  by 

</»,/»>      =  £     -T\f     /.(x)/*W<fc.  (6-14) 

7eF(S0)  «(7)  J**-* 

where  F(S°)  is  the  set  of  (triangular  or  quadrilateral,  depending  on  the  subdivision 
scheme)  faces  in  the  control  mesh  S°,  7  is  a  face  in  the  set  F(S°),  0(7)  is  the  area  of 
the  face  7  and  dx  is  the  usual  differential  area  in  St3. 
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su       J  S' 

Figure  6.5.  Subdivision  refinement  of  a  tetrahedron 


6.2.3     Biorthogonal  Surface  Wavelets  on  Arbitrary  Manifold 

The  biorthogonal  wavelets  discussed  here  can  be  developed  for  any  subdivision 
scheme  in  general.  However,  some  minor  modifications  are  needed  in  some  situations. 
The  construction  is  based  on  the  lifting  scheme  [93]  where  first  "lazy"  wavelets  at 
level  j  —  1  are  constructed  using  some  scaling  functions  at  level  j,  and  then  these 
wavelets  are  "lifted"  so  that  they  become  as  orthogonal  as  desired  to  the  scaling 
functions  at  level  j  -  1.  These  concepts  are  described  below  using  a  simple  example. 
A  more  detailed  discussion  on  this  construction  can  be  found  in  Stollnitz  et  al.  [92]. 

Let  3°,  as  shown  in  Fig.6.5,  be  a  (tetrahedral)  control  mesh  which  when  refined 
using  some  subdivision  rules  leads  to  a  smooth  surface  in  the  limit.  The  subdivision 
rules  refine  each  triangular  face  into  four  triangular  faces  by  introducing  new  vertices 
corresponding  to  each  edge  in  the  mesh  as  shown  in  Fig.6.5.  Let  the  smooth  surface 
in  the  limit  be  s  (denoted  by  S°°  in  Fig.6.5).  Let  p°  be  the  collection  of  four  control 
vertex  positions  defining  the  control  mesh  S°.  Now,  there  exists  scaling  functions 
0j(x),  x  G  S°,  associated  with  each  vertex  p,,  j  =  1, . . .  ,4,  such  that  they  span  the 
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space  V°(S°).  These  scaling  functions  depend  on  the  subdivision  scheme  involved 
and  may  or  may  not  have  analytic  expressions.  These  scaling  functions  are  collected 
in  J°(x)  =  [0?(x),  </>°(x),  0°(x),  0$(x)].  The  smooth  limit  surface  s  can  be  written 
as 

s  =  J°(x)  p°,  (6.15) 

where  x  e  5°.  The  control  mesh  Sl  is  obtained  by  subdividing  the  control  mesh  S" 
once.  The  ten  control  vertices  of  the  mesh  S1  can  be  categorized  into  two  classes  : 
(1)  four  control  vertices  corresponding  to  the  four  control  vertices  in  original  mesh 
S°  (the  filled  dots  in  Fig.6.5)  and  (2)  six  control  vertices  corresponding  to  the  six 
edges  in  original  mesh  S°  (the  unfilled  dots  in  Fig.6.5).  If  the  control  vertex  positions 
of  the  mesh  S1  are  collected  in  p1  such  that  the  "old"  vertices  corresponding  to  the 
vertices  in  the  original  mesh  precede  the  "new"  vertices  corresponding  to  the  edges 
in  the  original  mesh,  then 


py 


p1 


(6.16) 


where  P1  is  the  subdivision  matrix,  and  PlM  and  Pj,eu)  are  portions  of  the  subdivision 
matrix  that  encode  rules  to  obtain  the  "old"  and  "new"  vertices  respectively.  Basis 
functions  corresponding  to  the  vertices  in  p1  can  be  collected  in  Jl(x)  =  [3\id{x)  \ 
^liewi1)}  ~  [<Ai(x)>  •  •  ■  1 04(x)  |  ^(x), . . . ,  0}o(x)],  where  the  first  four  basis  func- 
tions are  associated  with  old  vertices  and  the  last  six  basis  functions  are  associated 
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with  new  vertices.  The  same  smooth  limit  surface  s  can  be  written  as 


s  =  J1(x)pl. 


(6.17) 


Evidently,  the  dimension  of  V°(Sa)  and  Vl(S°)  are  4  and  10  respectively.  In  the 
lazy  wavelet  construction,  a  complement  space  W°  (S°)  of  dimension  6  is  constructed 
with  the  set  of  basis  functions  corresponding  to  the  new  vertices  in  S1,  i.e.,  the 
wavelets  spanning  W°(S°)  are  W°(x)  =  [Vf  (x), . . . ,  <(x)]  =  [$to.  •  •  ■ .  <#o(x)]- 
This  construction  is  called  "lazy"  because  no  extra  computation  need  to  be  done  to 
derive  the  wavelets.  The  refinement  relations  for  the  lazy  wavelet  construction  can 
be  compactly  written  as 


[j°(x)  |  w°(x)]  =  py*)  I  JLJx)]  ppjL_  i  QL„ 


(6.18) 


where  P,1      =  P1  and 


[Plazy  I  Qlozj] 


0  and  I  denoting  zero  and  identity  matrix  respectively.  After  developing  the  synthesis 
filters  for  the  lazy  wavelets,  the  analysis  filters  can  be  derived  using  the  inverse  relation 
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and  are  given  by 

r  a1    " 

P1 

0 

-i 

(p^r1 

0 

."«•»». 

P1 

*  new 

I 

pi     /pi    y 

nem\x   o/d/ 

1    I 

(6.19) 


The  lazy  wavelet  construction  developed  here  can  be  generalized  for  any  subdivision 
level  j  in  a  straight-forward  manner.  The  synthesis  and  analysis  filters  at  level  j  for 
a  subdivision  scheme  can  be  written  as 


\?lazy  I  Qiasy] 


PL    o 


(6.20) 


and 


lazy 


BL, 


(P50„ 


(6.21) 


respectively,  where  P£H  encodes  subdivision  rules  for  obtaining  "old"  vertices  at  j- 
th  level  corresponding  to  the  vertices  in  the  {j  -  l)-th  level  control  mesh  and  Pieul 
encodes  rules  for  obtaining  the  other  "new"  vertices  in  the  j'-th  level  control  mesh. 

It  may  be  noted  that  the  lazy  wavelet  construction  leads  to  sparse  synthesis 
matrices.  However,  the  corresponding  analysis  filters  may  not  be  sparse  in  general  as 
inverse  of  P^((i  is  not  necessarily  sparse.  This  implies  that  even  though  the  synthesis 
can  be  done  in  linear  time  with  respect  to  the  number  of  vertices,  the  complexity  of 
the  analysis  depends  on  how  fast  the  inverse  can  be  computed.  The  synthesis  can  not 
be  done  in  linear  time  in  case  of  approximating  subdivision  schemes.  However,  PJold 
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Wavelet  associated  with  the  edge 
connecting  a  and  b 

Lazy  :  the  scaling  function  associated  with 
the  introduced  vertex  at  higher  level 

Lifted  :  the  scaling  function    associated  with 
the  introduced  vertex  at  higher  level 
-  (  weighted  combination  of  the  scaling 
functions  associated  with  the  solid 
vertices) 

Figure  6.6.  Wavelet  construction. 


is  an  identity  matrix  for  interpolatory  subdivision  schemes,  and  hence  both  synthesis 
and  analysis  can  be  done  in  linear  time  for  interpolatory  subdivision  schemes.  The 
simplified  expressions  of  the  synthesis  and  analysis  filters  in  case  of  interpolatory 
subdivision  scheme  are  given  by 


l^lazy  1  Qlazy]     ~ 

I           0 

pj      i 

new 

and 

"■lazy 

I         0 

.Mr. 

-Pj        I 

new 

respectively. 

(6.22) 


(6.23) 
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The  lazy  wavelets  are  easy  to  construct.  However,  a  lower  resolution  mesh 
obtained  using  the  analysis  filters  of  the  lazy  wavelet  construction  is  far  from  the 
least-squares-best  approximation.  Lazy  wavelets  at  some  level  (j  —  1)  are  nowhere  or- 
thogonal to  V'~l.  The  lifting  scheme  [93]  improves  the  lazy  wavelets  at  level  (j'-l)  by 
subtracting  a  liner  combination  of  the  coarser  level  scaling  functions  which  overlap  to 
a  given  wavelet.  The  resulting  wavelet  function  becomes  somewhat  more  orthogonal 
to  the  coarser  level  scaling  functions  than  the  original  lazy  wavelet.  Mathematically, 
a  "lifted"  wavelet  at  level  (j  —  1)  can  be  expressed  as 


I 
=    ^L-»W-S*W-1W.  (6-24) 

i 


where  i-th  lazy  wavelet  function  at  level  (j  —  1)  is  the  scaling  function  associated 
with  the  i-th  "new"  vertex  at  level  j,  and  /  is  restricted  to  a  few  vertices  at  level 
(j  —  1)  in  the  neighborhood  of  the  i-th  new  vertex  at  level  j  (Fig.6.6).  In  order  to  get 
the  values  of  sjt  for  the  lifted  wavelet  t,  a  linear  system  of  equations  is  derived  and 
solved  in  a  least  squares  sense  by  imposing  the  constraints  <  ^^(x),  0J_1(x)  >  = 
0  for  all  k,  such  that  the  supports  of  ^i,/(ij(x)  and  </^~  (x)  overlap.  The  number  of 
non-zero  entries  is  in  sj  i  needs  to  be  controlled.  A  large  number  of  non-zero  entries 
will  lead  to  a  lifted  wavelet  function  that  is  nearly  orthogonal  to  the  coarser  level 
scaling  functions,  at  the  same  time  the  multiresolution  analysis  will  become  more 
time  consuming,  and  hence  a  tradeoff  needs  to  be  decided  as  far  as  orthogonality 
and  time  of  computation  is  concerned.  In  Stollnitz  et  al.  [92],  "k-disk  wavelets"  are 
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constructed  for  triangular  mesh-based  subdivision  schemes  where  the  non-zero  entries 
of  sjj  are  restricted  to  the  vertices  at  level  (j'-l)  which  are  within  the  fc-neighborhood 
of  the  i-th  "new"  vertex  at  level  j.  The  reader  is  referred  to  Stollnitz  et  al.  [92]  for 
the  details  of  the  construction.  Once  the  non-zero  values  of  sf  {  are  computed,  they 
can  be  assembled  to  a  matrix  S> .  Then,  the  expressions  for  the  synthesis  and  analysis 
filters  of  the  lifted  wavelet  scheme  can  be  written  as 


[P'lift  I  Qliftl     -     [P/ass,  I  Qfazy        PtmyS'l 


(6.25) 


and 


Auft 


'lift 


ALy  +  SJBL„ 


alazy 


(6.26) 


respectively. 

6.3     Multiresolution  Representation 

The  multiresolution  analysis  on  arbitrary  manifold  was  developed  to  obtain  good 
low  resolution  approximations  of  complicated  meshes.  This  is  a  top-down  technique 
where  one  starts  with  a  complicated  arbitrary  topology  mesh  and  goes  on  obtain- 
ing lower  resolution  meshes  using  analysis  filters.  The  complicated  mesh  on  which 
multiresolution  analysis  is  applied  can  be  reconstructed  using  the  synthesis  filters. 
In  this  dissertation,  the  multiresolution  technique  is  used  in  a  novel  bottom-up  ap- 
proach where  an  evolving  control  mesh  of  a  dynamic  subdivision  surface  is  subdivided 
at  equilibrium  to  obtain  more  degrees  of  freedom,  and  then,  the  resulting  control  mesh 
is  treated  as  the  previous  control  mesh  and  a  collection  of  detail  (wavelet)  coefficients. 
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Of  course,  the  detail  coefficients  are  zero  when  an  control  mesh  is  replaced  by  another 
control  mesh  obtained  by  a  pure  subdivision  step  on  the  previous  mesh,  but  becomes 
non-zero  as  the  new  control  mesh  (represented  as  previous  control  mesh  and  wavelet 
coefficients)  evolves  over  time  in  response  with  the  synthesized  force  application.  This 
idea  is  further  illustrated  in  the  rest  of  this  section. 

Let  s  be  a  smooth  limit  surface  obtained  via  infinite  number  of  subdivision 
steps  on  an  initial  mesh  S°,  whose  vertex  positions  are  collected  in  p°.  As  mentioned 
earlier,  the  limit  surface  can  be  written  as 

s  =  J°p°,  (6.27) 

where  J°  is  a  collection  of  basis  functions.  If  the  control  mesh  S1  is  obtained  by 
subdividing  S°  once,  then  the  same  limit  surface  can  be  expressed  as 

s  =  jy,  (6.28) 

where  p1  is  the  collection  of  vertex  positions  of  the  mesh  S1  and  J1  is  the  collection 
of  corresponding  basis  functions.  Since  p1  =  P'p0,  the  limit  surface  s  can  also  be 
written  as 

s  =  J'Py.  (6.29) 

The  difference  between  Eqn.6.27  and  Eqn.6.28  is  that  the  limit  surface  s  has  more 
degrees  of  freedom  (control  vertices)  for  representation  in  the  latter  in  comparison 
with  the  former.   However,  in  Eqn.6.29,  even  though  the  limit  surface  is  expressed 
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using  the  same  basis  functions  as  that  of  S1,  the  degrees  of  freedom  remains  the  same 
as  that  of  S°.  Eqn.6.29  can  be  modified  in  the  following  manner  so  that  it  uses  same 
basis  functions  and  same  number  of  degrees  of  freedom  as  that  of  S  . 


s  =  J1  [P'Q 


(6.30) 


where  q°  is  the  collection  of  wavelet  coefficients  whose  values  are  zero,  P1  and  Q1 
are  the  synthesis  matrices  P,1^,  andQji/(  respectively.  The  above  formulation  yields 
a  multiresolution  representation  of  the  control  mesh  S1  obtained  via  one  subdivision 
step  on  the  previous  control  mesh  S°.  The  mesh  S1  is  represented  as  a  mesh  at  level 
0  along  with  the  detail  (wavelet)  coefficients  at  level  0.  As  mentioned  earlier,  these 
coefficients  would  become  non-zero  when  the  smooth  limit  surface  deforms  over  time 
due  to  synthesized  force  application.  Similarly  after  one  more  subdivision  step,  the 
new  control  mesh  S2  can  be  viewed  as  control  vertex  positions  at  level  0,  wavelet 
coefficients  at  level  0  and  wavelet  coefficients  at  level  1,  and  can  be  written  as 


s    =    J    p 


J2  [P2Q2] 


J2  [P2Q2 
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P1 
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q° 
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q' 

(6.31) 
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where  J2  is  the  collection  of  basis  functions  at  level  2,  p2,  p1  and  p°  are  the  control 
vertex  positions  at  level  2,  1  and  0  respectively,  q1  and  q°  are  the  wavelet  coefficients 
at  level  1  and  0  respectively,  I  is  the  identity  matrix  and  0  is  the  zero  matrix.  This 
idea  can  be  generalized  for  a  control  mesh  SJ  obtained  after  j'-th  subdivision  and  the 
corresponding  expression  can  be  written  as 


=     TJ    AJ 


33  A1  pi, 


(6.32) 


where  (pj)T  =  [p°q° . . .  qJ'~  V_1]>  and 
A'  =  [PJQ 


pj-i    Qj-i    o 
0  0       I 


pj-2      QJ-2      0      0 

0  Oil 


P°    Q°    0 
0      0      1 


The  multiresolution  dynamics  is  developed  using  this  formulation  in  the  next 
section. 

6.4     Dynamics 

The  smooth  limit  surface  s  can  be  made  dynamic  if  the  control  vertex  positions 
at  level  0  and  wavelet  coefficients  at  level  1, 2, ...,  (j  -  1)  are  functions  of  time.  The 
velocity  of  the  smooth  surface,  controlled  by  the  mesh  Sj  obtained  through  j  steps 
of  model  subdivision,  is  given  by 


i(x,p>)  =  J'(x)AVr. 


(6.33) 
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where  the  overstruck  dot  denotes  the  time  derivative  and  x  e  S*.  The  Lagrangian 
motion  equation  (refer  to  Eqn.2.5)  can  be  derived  in  a  similar  fashion  as  mentioned 
in  earlier  chapters,  and  is  given  by 

M*#  +  Djp]T  +  KJp?  =  f/,  (6.34) 

where  MJ,  DJ  and  KJ  are  the  j-th  level  mass,  damping  and  stiffness  matrices  respec- 
tively and  fj  is  the  generalized  force  vector  at  level  j.  If  /i(x)  is  the  mass  density  of 
the  subdivision  surface  model,  then  the  mass  matrix  at  j-th  level  is  given  by 


Mj  =  f      /i(x)(AJ)r(JJ'(x))TJJ(x)AJrfx. 
Jxes> 


The  damping  matrix  EV  can  be  derived  in  a  similar  fashion,  while  the  stiffness  ma- 
trix K^  can  be  obtained  following  the  techniques  described  in  Section  4.3.2.  The 
generalized  force  vector  at  level  j  can  be  written  as 


Vv=l       (AJ)V(x))Tf(x,0<*x. 


6,5     Implementation  Details 

Multiresolution  representation  of  the  evolving  control  mesh  achieves  coarse-to- 
fine  as  well  as  fine-to-coarse  manipulation  of  the  smooth  limit  surface.  For  example, 
the  user  starts  with  a  deformable  smooth  surface  with  a  simple  control  mesh  5° 
and  directly  manipulates  the  limit  surface  by  applying  synthesized  forces.  When  an 


130 


equilibrium  is  obtained,  the  user  can  increase  the  degrees  of  freedom  by  switching  to 
a  control  mesh  Sl  which  is  theoretically  obtained  by  applying  one  subdivision  step 
on  S°,  but  represented  as  a  collection  of  vertex  positions  in  S°  along  with  wavelet 
coefficients  at  level  0  needed  to  reconstruct  51.  The  wavelet  coefficients  are  zero  at 
the  point  of  switching,  but  becomes  non-zero  as  the  control  mesh  S1,  represented 
as  mentioned  above,  evolves  over  time.  The  user  can  further  increase  the  degrees 
of  freedom  to  obtain  more  localized  effect  by  introducing  control  mesh  S2  which  is 
represented  as  control  vertex  positions  of  mesh  5°,  along  with  wavelet  coefficients 
at  level  0  and  1.  Thus  the  user  can  have  coarse-to-fine  manipulation  of  the  smooth 
limit  surface.  This  facility  is  also  present  in  the  multilevel  dynamics  approach  as 
mentioned  earlier.  However,  the  multilevel  dynamics  approach  does  not  support  fine- 
to-coarse  manipulation  -  it  fails  if  the  user  wants  to  apply  synthesized  forces  on  a 
coarser  control  mesh  at  level  j  after  moving  up  to  level  j  +  k.  This  can  be  achieved 
in  multiresolution  approach  by  simply  using  the  dynamic  equation  of  motion  at  level 
j,  and  applying  synthesized  forces  on  a  smooth  limit  surface  which  is  obtained  by 
multiresolution  synthesis  using  "time  varying"  control  vertex  positions  at  level  0 
and  wavelet  coefficients  at  level  0, 1, . . . ,  (j  —  1)  and  "static"  wavelet  coefficients  at 
level  j,  (j  +  1), . . . ,  {j  +  k  —  1).  If  the  user  switches  back  to  level  (j  +  k),  then 
the  control  vertex  positions  at  level  0  as  well  as  all  the  wavelet  coefficients  at  level 
0, 1, . . . ,  (j  -  1),  j,  (j  +  1), . . . ,  (j  +  k  —  1)  become  function  of  time,  and  the  system 
switches  back  to  the  motion  equation  of  level  (j  +  k). 
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The  actual  implementation  differs  slightly  from  the  formulation  to  achieve  effi- 
ciency in  coarse-to-fine  and  fine-to-coarse  manipulation.  In  the  implemented  system, 
the  user  starts  out  with  a  smooth  limit  surface  which  has  a  simple  control  mesh  with 
very  few  degrees  of  freedom.  The  model  grows  and  when  the  equilibrium  is  achieved, 
one  step  of  subdivision  yields  larger  degrees  of  freedom.  This  larger  degrees  of  free- 
dom are  the  control  vertex  positions  in  this  finer  resolution  and  not  control  vertex 
positions  in  lower  resolution  along  with  wavelet  coefficients  in  that  lower  resolution. 
The  model  can  grow  further,  and  if  more  subdivision  step  is  needed  to  increase  the 
degrees  of  freedom,  it  is  handled  in  a  similar  fashion.  If  at  any  certain  point  of  time 
the  user  needs  to  go  back  to  a  coarser  level  control  mesh,  wavelet  decomposition  is 
done  on  the  higher  resolution  base  mesh.  This  decomposition  yields  vertex  positions 
in  the  coarser  level  mesh  (which  becomes  the  new  base  mesh)  and  non-zero  wavelet 
coefficients.  The  lower  resolution  mesh  evolves  with  time  due  to  the  force  applied  on 
the  limit  surface.  It  may  be  noted  that  the  limit  surface  in  this  scenario  is  obtained  by 
doing  wavelet  reconstruction  with  non-zero  wavelet  coefficients  at  some  levels.  The 
user  can  choose  further  coarser  level  control  mesh  by  repeating  the  same  process. 

In  the  implemented  system,  the  user  can  achieve  global  deformation  on  a  de- 
tailed mesh  by  going  down  couple  of  levels  using  wavelet  decomposition,  retaining 
the  wavelet  coefficients,  and  applying  force  on  the  limit  surface  (obtained  through 
wavelet  reconstruction  process)  and  letting  the  coarser  level  mesh  evolve.  Similarly, 
to  achieve  local  deformation  on  a  coarser  mesh,  the  user  can  move  up  couple  of  steps 
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by  doing  subdivision  (and  adding  non-zero  wavelet  coefficients  at  higher  levels,  if 
present)  and  applying  synthesized  forces  in  the  region  of  interest. 

The  formulation  using  wavelet  representation  has  mathematical  niceties,  but 
inefficient  for  implementation  purposes.  The  synthesized  force  is  applied  on  the  limit 
surface.  To  evaluate  the  limit  surface,  wavelet  reconstruction  needs  to  be  done  at  each 
time  step  starting  from  the  original  control  mesh  upto  a  specified  level  if  an  explicit 
wavelet  representation  of  the  evolving  control  mesh  is  maintained  as  mentioned  in 
the  formulation.  On  the  other  hand,  if  the  evolving  mesh  is  control  vertex  positions 
at  that  resolution,  then  the  wavelet  reconstruction  starts  from  that  resolution  for 
evaluating  the  limit  surface.  Wavelet  decomposition  needs  to  be  done  while  going 
down  from  finer  to  coarser  resolution,  but  this  happens  only  once  in  a  while,  and  not 
at  each  time  step.  Therefore,  the  implemented  version  is  much  more  efficient.  It  may 
be  noted  that  the  motion  equation  at  any  level  is  implemented  using  finite  element 
techniques  discussed  in  earlier  chapters  for  various  subdivision  schemes. 


CHAPTER  7 
SYSTEM  ARCHITECTURE 


The  architecture  of  the  prototype  system  is  described  in  this  chapter.  The  basic 
building  blocks  of  the  system  are  different  modules  performing  different  tasks.  The 
overall  architecture  of  the  system  is  shown  in  Fig.7.1. 

The  user  first  provides  the  input  about  the  geometry  and  topology  of  the  initial 
mesh.  Typically  these  information  consist  of  some  type  of  connectivity  information 
and  control  vertex  positions  in  3D.  The  topological  information  processing  module 
derives  other  necessary  topological  information  from  the  given  input  data,  and  passes 
all  the  information  to  the  subdivision  module.  The  subdivision  module  performs  pre- 
specified  number  of  subdivision  steps  on  the  initial  mesh,  and  then  the  finite  element 
analysis  module  computes  different  elemental  matrices  needed  by  the  update  engine. 
The  user  can  apply  different  types  of  synthesized  forces  on  the  limit  surface,  and 
the  model  is  constantly  updated  to  obtain  an  equilibrium  position.  The  updated 
information  is  passed  to  the  display  module  which  renders  the  deforming  smooth 
surface  on  the  screen.  The  user  is  also  allowed  to  choose  different  model  parameter 
values  like  mass  density,  damping  density,  spring  constants  etc.  Details  of  the  specific 
modules  are  provided  in  the  rest  of  this  chapter. 
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Figure  7.1.  System  architecture. 
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7.1     Topological  Information  Processing  Module 

In  the  implemented  system,  the  user  needs  to  provide  the  connectivity  of  the 
initial  mesh  by  listing  the  vertices  of  every  face.  The  vertex  information  of  a  face 
must  be  provided  in  an  orderly  fashion  (either  clockwise  or  anti-clockwise).  The 
topological  information  processing  module  derives  other  necessary  information  like 
vertex  connectivity  information  of  the  edges,  edge  connectivity  information  of  the 
faces  etc.  from  the  input  provided  by  the  user. 

7.2     Subdivision  Module 

This  module  applies  different  subdivision  rules  to  refine  an  input  control  mesh. 
Currently,  Catmull-Clark  and  modified  butterfly  subdivision  schemes  are  supported 
in  the  system.  Other  subdivision  schemes  can  also  be  easily  incorporated.  The 
subdivision  scheme  for  a  specific  application  is  chosen  by  the  user.  However,  there 
are  certain  inherent  restrictions  depending  on  the  chosen  subdivision  scheme.  For 
example,  butterfly  subdivision  scheme  can  be  chosen  only  if  the  initial  (control) 
mesh  has  triangular  faces.  Once  the  user  provides  input  about  the  initial  mesh 
and  the  chosen  subdivision  scheme,  this  module  then  subdivides  the  initial  mesh 
a  pre-specified  number  of  times  to  obtain  a  high  resolution  mesh  which  effectively 
approximates  the  smooth  limit  surface. 

7.3     Finite  Element  Analysis  Module 

The  finite  element  analysis  module  determines  the  type  of  finite  element  to 
be  used  and  the  total  number  of  elements  describing  the  smooth  limit  surface.  It 
computes  all  elemental  matrices  such  as  mass,  damping  and  stiffness  matrices,  as 
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well  as  determines  the  set  of  vertices  in  the  initial  mesh  controlling  each  individual 
element.  This  module  also  provides  the  information  about  the  global  mass,  damping 
and  stiffness  matrices  without  actually  assembling  them  from  individual  elemental 
matrices. 

7.4     Force  Synthesis  Module 

The  user  can  apply  various  types  of  forces  on  the  smooth  limit  surface.  Currently 
spring  forces,  balloon  forces,  and  image  gradient-based  forces  are  being  supported. 
Other  types  of  forces  can  also  be  easily  incorporated.  The  user  applies  spring  forces 
on  the  limit  surface  by  providing  points  in  3D  from  which  the  springs  need  to  be 
attached.  The  3D  point  specification  is  done  either  interactively  from  mouse  input  or 
from  files.  To  apply  balloon  forces,  the  user  needs  to  specify  whether  the  ballooning 
force  is  to  inflate  the  model  or  to  deflate  the  model.  Volume  images  are  provided 
as  input  to  apply  image  gradient-based  forces.  A  detailed  discussion  on  these  force 
application  methods  was  provided  in  Section  3.3.4. 

Lh Update  Engine 

After  the  user  provides  all  the  necessary  information  about  the  model  and  the 
force  application  method,  the  update  engine  constantly  updates  the  dynamic  sub- 
division surface  model  to  achieve  an  equilibrium  position.  The  discretized  second 
order  differential  equation  of  motion  leads  to  solving  a  large  sparse  linear  system 
of  equations.  Conjugate  gradient  [34,  75],  an  iterative  linear  system  solver,  is  used 
for  this  purpose.   The  user  can  choose  either  multilevel  dynamics  or  multiresolution 
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dynamics  approach,  and  the  update  engine  accordingly  chooses  algorithms  when  the 
degrees  of  freedom  of  the  evolving  surface  model  is  changed. 

7.6    Display  Module 

The  update  engine  constantly  updates  the  dynamic  model,  and  the  relevant 
geometric  and  topological  information  is  passed  to  the  display  module  for  rendering 
the  evolving  smooth  surface.  The  display  module  uses  OpenGL  graphics  libraries 
for  rendering  purposes.  It  also  serves  as  the  input  module  for  mouse-based  user 
interactions. 


CHAPTER  8 
APPLICATIONS 


The  proposed  dynamic  framework  extends  the  application  areas  of  subdivision 
surfaces  beyond  the  traditional  modeling  domain.  The  framework  enhances  the  suit- 
ability of  subdivision  surfaces  for  geometric  modeling.  At  the  same  time,  the  proposed 
dynamic  models  can  be  used  for  recovery  of  shapes  from  range  and  volume  data,  non- 
rigid  motion  estimation,  multiresolution  editing  and  natural  terrain  modeling.  The 
applicability  of  the  dynamic  model  in  these  areas  is  further  illustrated  in  the  rest  of 
this  chapter. 

8.1     Geometric  Modeling 

The  proposed  dynamic  framework  for  subdivision  surfaces  can  be  used  to  rep- 
resent a  wide  variety  of  arbitrary  topological  shapes  with  remarkable  ease.  In  a 
typical  geometric  modeling  application  using  dynamic  subdivision  surfaces,  the  user 
can  specify  any  mesh  as  the  initial  (control)  mesh,  and  the  corresponding  limit  surface 
can  be  sculpted  interactively  by  applying  synthesized  forces.  Any  dynamic  subdivi- 
sion surface  model  can  be  used  for  the  modeling  purpose.  The  examples  shown  in 
this  section  use  dynamic  Catmull-Clark  and  dynamic  modified  butterfly  subdivision 
surface  models. 
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Figure  8.1.  (a),  (b),  (c)  and  (d)  :  Initial  shapes  (obtained  applying  Catmull-Clark 
subdivision  rules  on  control  meshes);  (e),  (f),  (g)  and  (h)  :  the  corresponding  modified 
shapes  after  interactive  force  application. 
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First,  dynamic  Catmull-Clark  surface  models  as  presented  in  Chapter  3  are  used 
for  shape  modeling.  The  special  elements  in  the  limit  surface  are  darkly  shaded  and 
the  normal  elements  are  lightly  shaded  in  all  the  modeling  examples.  It  may  be 
noted  that  the  formulation  presented  in  Chapter  5,  which  will  have  only  one  type 
of  finite  element  in  the  limit  surface,  can  be  used  as  well.  In  Fig.8.1,  several  initial 
surfaces  obtained  from  different  control  meshes  using  Catmull-Clark  subdivision  rules 
and  the  corresponding  modified  surfaces  after  interactive  spring  force  application  are 
shown.  To  change  the  shape  of  an  initial  surface,  springs  are  attached  from  different 
points  in  3D  to  the  nearest  points  on  the  limit  surface  such  that  the  limit  surface 
deforms  towards  these  points  generating  the  desired  shape.  In  Fig.8.1  (a),  an  open 
surface  defined  by  an  initial  mesh  of  61  vertices  and  45  faces  is  shown.  The  mesh 
has  one  extraordinary  vertex  of  degree  5.  This  limit  surface  is  modified  by  applying 
spring  forces  interactively,  and  the  modified  surface  is  depicted  in  Fig.8.1  (e).  A 
torus,  defined  by  an  initial  mesh  of  32  vertices  and  32  faces,  and  its  modified  shape  is 
shown  in  Fig.8.1(b)  and  (f)  respectively.  The  initial  mesh  of  the  smooth  limit  surface 
shown  in  Fig.8.1(c)  has  544  faces  and  542  vertices,  8  of  them  are  extraordinary 
vertices  of  degree  5.  The  limit  surface  is  modified  interactively  by  applying  spring 
forces  from  various  points  in  3D  and  the  modified  shape  is  depicted  in  Fig.8.1(g). 
Note  that  the  extent  of  deformation  has  been  interactively  controlled  by  varying  the 
stiffness  of  the  attached  springs.  The  upper  portion  of  the  limit  surface  has  been 
deformed  by  applying  spring  forces  of  higher  magnitude,  whereas  the  lower  portion 
has  been  modified  by  applying  spring  forces  of  lower  magnitude.  The  spread  of  the 
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Figure  8.2.  (a),  (b),  (c)  and  (d)  :  Initial  shapes  (obtained  applying  modified  butterfly 
subdivision  rules  on  control  meshes);  (e),  (f),  (g)  and  (h)  :  the  corresponding  modified 
shapes  after  interactive  sculpting  via  force  application. 


deformation  effect  is  clearly  larger  in  the  latter  case  for  obvious  reasons.  Finally,  a 
flat  sheet  defined  by  an  initial  mesh  of  64  faces  and  81  vertices,  shown  in  Fig.8.1(d), 
is  deformed  interactively  to  obtain  the  hat-like  shape  shown  in  Fig.8.1(h). 

Next,  shape  modeling  examples  using  the  dynamic  modified  butterfly  subdivi- 
sion surface  model  are  shown.  The  limit  surface  here  consists  of  a  single  type  of 
smooth  triangular  finite  element  patches,  irrespective  of  the  number  of  extraordinary 
vertices  in  the  control  mesh.  In  Fig.8.2,  several  initial  surfaces  obtained  from  different 
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control  meshes  using  the  modified  butterfly  subdivision  rules  and  the  corresponding 
modified  surfaces  after  interactive  sculpting  are  shown.  As  with  the  dynamic  Catmull- 
Clark  subdivision  surfaces,  the  user  can  attach  springs  from  different  points  in  3D  to 
the  nearest  points  on  the  limit  surface  and  the  limit  surface  deforms  towards  these 
points  to  generate  the  desired  shape.  The  initial  mesh  of  the  smooth  surface  shown 
in  Fig.8.2(a)  has  125  faces  and  76  vertices  (degrees  of  freedom),  which  is  deformed 
to  the  smooth  shape  shown  in  Fig.8.2(e)  by  interactive  spring  force  application.  The 
initial  mesh  of  the  closed  cube-like  shape  in  Fig.8.2(b)  has  24  faces  and  14  vertices. 
This  cube-like  shape  is  deformed  to  the  shape  shown  in  Fig.8.2(f).  The  one  hole  torus 
in  Fig.8.2(c)  and  the  corresponding  modified  shape  in  Fig.8.2(g)  have  initial  meshes 
with  64  faces  and  32  vertices.  A  two  hole  torus  with  a  control  mesh  of  272  faces 
and  134  vertices,  shown  in  Fig.8.2(d),  is  dynamically  sculpted  to  the  shape  shown  in 
Fig.8.2(h). 

&2 Shape  Recovery  from  Range  and  Volume  Data 

The  dynamic  subdivision  surface  models  can  recover  the  underlying  shape  of 
a  given  range  or  volume  data  set  in  a  hierarchical  fashion.  The  initialized  model 
deforms  under  the  influence  of  the  synthesized  forces  from  the  range  or  volume  data. 
When  an  approximate  shape  is  recovered,  a  new  control  mesh  can  be  obtained  by 
doing  one  step  of  subdivision  on  the  initial  mesh  thereby  increasing  the  degrees  of 
freedom  to  represent  the  same  limit  surface,  and  a  better  fit  to  the  given  range  data 
set  can  be  achieved.  It  may  be  noted  that  the  model  initialization  is  interactive, 
and  the  initialized  model  can  have  any  control  mesh  of  the  desired  genus.  However, 
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Figure  8.3.  (a)  Range  data  of  a  bulb  along  with  the  initialized  model,  (b)  an  in- 
termediate stage  of  evolution,  and  (c)  the  fitted  dynamic  Catmull-Clark  subdivision 
surface  model. 


an  initial  mesh  with  few  degrees  of  freedom  usually  performs  better  in  terms  of  the 
compact  representation  of  the  underlying  shape. 

Any  dynamic  subdivision  surface  model  can  be  used  for  shape  recovery.  Next, 
three  examples  of  shape  recovery  from  range  data  sets  are  shown.  The  first  two 
examples  use  dynamic  Catmull-Clark  subdivision  surface  model  and  the  third  one 
uses  dynamic  modified  butterfly  subdivision  surface  model.  The  error  in  fit,  which 
is  defined  as  the  maximum  distance  between  a  data  point  and  the  nearest  point  on 
the  limit  surface  as  a  percentage  of  the  diameter  of  the  smallest  sphere  enclosing  the 
object,  is  approximately  3%  in  all  the  experiments  with  range  data.  The  time  of 
dynamic  evolution  for  fitting  the  range  data  sets  used  in  the  experiments  is  approxi- 
mately 3  minutes  in  a  SGI  02  workstation. 
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Figure  8.4.  (a)  Range  data  of  an  anvil  along  with  the  initialized  model,  (b)  an 
intermediate  stage  of  evolution,  and  (c)  the  fitted  dynamic  Catmull-Clark  subdivision 
surface  model. 


(b) 


to 


Figure  8.5.  (a)  Range  data  of  a  head  along  with  the  initialized  model,  (d)  the  fitted 
dynamic  butterfly  subdivision  model,  and  (c)  visualization  of  the  shape  from  another 
view  point. 
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The  examples  shown  in  Fig.8.3  and  Fig.8.4  involve  dynamic  Catmull-Clark  sub- 
division surface.  The  same  shading  convention  for  normal  and  special  elements  as 
mentioned  earlier  is  used.  In  both  of  these  examples,  the  initialized  model  had  96 
faces  and  98  vertices,  8  of  them  being  extraordinary  vertices  of  degree  3.  The  final 
fitted  model,  obtained  through  one  step  of  subdivision,  has  a  control  polygon  of  384 
faces  with  386  vertices.  In  Fig.8.3,  the  dynamic  Catmull-Clark  subdivision  surface 
model  is  fitted  to  laser  range  data  acquired  from  multiple  views  of  a  light  bulb.  Prior 
to  applying  the  model  fitting  algorithm,  the  data  were  transformed  into  a  single  refer- 
ence coordinate  system.  The  model  was  initialized  inside  the  1000  range  data  points 
on  the  surface  of  the  bulb.  In  the  next  experiment,  the  dynamic  Catmull-Clark  sub- 
division surface  model  is  fitted  to  an  anvil  data  set  (Fig.8.4).  The  data  set  has  2031 
data  points.  It  may  be  noted  that  the  final  shape  with  3%  error  tolerance  uses  very 
few  control  points  for  representation  in  comparison  with  the  number  of  data  points 
present  in  the  original  range  data  set. 

In  the  last  example  with  range  data  set,  the  shape  of  a  human  head  is  recovered 
from  a  range  data  set  using  dynamic  modified  butterfly  subdivision  surface  model 
(Fig.8.5).  The  initialized  model  has  a  control  mesh  comprising  of  24  triangular  faces 
and  14  vertices  whereas  the  control  mesh  of  the  fitted  model  has  384  triangular 
faces  and  194  vertices.  The  range  data  set  has  1779  points  in  3D,  and  a  compact 
representation  of  the  recovered  shape  using  few  degrees  of  freedom  in  comparison 
with  original  data  set  is  obtained  in  this  example  as  well. 
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Figure  8.6.  (a)  A  slice  from  a  brain  MM,  (b)  initialized  model  inside  the  region  of 
interest  superimposed  on  the  slice,  (c)  the  fitted  model  superimposed  on  the  slice, 
and  (d)  a  3D  view  of  the  dynamic  Catmull-Clark  subdivision  surface  model  fitted  to 
the  cerebellum. 
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Figure  8.7.  (a)  Data  points  identifying  the  boundary  of  the  caudate  nucleus  on  a  MM 
slice  of  human  brain,  (b)  data  points  (from  all  slices)  in  3D  along  with  the  initialized 
model,  and  (c)  the  fitted  dynamic  butterfly  subdivision  model. 
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The  application  of  the  dynamic  subdivision  surface  models  to  anatomical  shape 
recovery  from  3D  volumetric  MRI  data  is  shown  in  the  next  two  examples.  The 
first  one  uses  dynamic  Catmull-Clark  model  and  the  next  one  uses  dynamic  modified 
butterfly  subdivision  surface  model. 

A  dynamic  Catmull-Clark  subdivision  surface  model  is  fitted  to  a  cerebellum 
(a  cortical  structure  in  brain)  given  an  input  of  30  sagittal  slices  from  a  MR  brain 
scan.  As  in  the  examples  with  range  data,  the  initialized  model  had  96  faces  and  98 
vertices,  8  of  them  being  extraordinary  vertices  of  degree  3.  The  final  fitted  model, 
obtained  through  one  step  of  subdivision,  has  a  control  polygon  of  384  faces  with  386 
vertices.  Fig.8.6(a)  depicts  a  slice  from  this  MRI  scan  and  the  model  initialization 
is  shown  in  Fig.8.6(b).  Continuous  image-based  forces  are  applied  to  the  model  and 
the  model  deforms  under  the  influence  of  these  forces  until  maximum  conformation 
to  the  boundaries  of  the  desired  cerebellum  shape.  The  final  fitted  model  is  shown 
in  Fig.8.6(c).  A  3D  view  of  the  fitted  model  is  depicted  in  Fig.8.6(d). 

In  the  next  example,  the  shape  extraction  of  a  caudate  nucleus  (another  cortical 
structure  in  human  brain)  is  presented  from  64  MRI  slices,  each  of  size  (256,256). 
Fig.8.7(a)  depicts  a  slice  from  this  MRI  scan  along  with  the  points  placed  by  an 
expert  neuroscientist  on  the  boundary  of  the  shape  of  interest.  Fig. 8. 7(b)  depicts  the 
interactively  placed  sparse  set  of  data  points  (placed  in  some  of  the  slices  depicting  the 
boundary  of  the  shape  of  interest)  in  3D  along  with  the  initialized  model.  Note  that 
points  had  to  be  interactively  placed  on  the  boundary  of  the  caudate  nucleus  in  MR 
slices  lacking  image  gradients  which  delineate  the  caudate  from  the  surrounding  tissue 
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in  the  image.  Continuous  image-based  forces  as  well  as  spring  forces  are  applied  to 
the  model  and  the  model  deforms  under  the  influence  of  these  forces  until  maximum 
conformation  to  the  boundaries  of  the  desired  caudate  shape.  The  final  fitted  model 
in  3D  is  shown  in  Fig. 8. 7(c).  The  initialized  dynamic  butterfly  subdivision  surface 
model  has  a  control  mesh  comprising  of  24  triangular  faces  and  14  vertices  whereas 
the  control  mesh  of  the  fitted  model  has  384  triangular  faces  and  194  vertices. 

8.3     Non-rigid  Motion  Estimation 

The  dynamic  subdivision  surface  models  can  be  used  effectively  to  estimate  the 
motion  of  a  non-rigid  object  from  a  given  time  sequence  of  range  or  volume  data.  An 
example  is  shown  in  Fig.8.8  where  the  motion  of  the  left-ventricular  chamber  of  a 
canine  heart  is  tracked  over  a  complete  cardiac  cycle  using  dynamic  modified  butterfly 
subdivision  surface  model.  The  data  set  comprised  of  16  3D  CT  images,  with  each 
volume  image  having  118  slices  of  128  x  128  pixels.  First,  the  shape  is  recovered  from 
one  data  set  using  image-based  (gradient)  as  well  as  point-based  forces.  As  before, 
the  initialized  model  has  a  control  mesh  comprising  of  24  triangular  faces  and  14 
vertices  whereas  the  control  mesh  of  the  fitted  model  has  384  triangular  faces  and 
194  vertices.  Once  the  shape  is  recovered  from  one  data  set,  this  fitted  model  is 
used  as  an  initialization  for  the  next  data  set  to  track  the  shape  of  interest.  The 
snapshots  from  motion  tracking  are  shown  in  Fig.8.8  for  the  16  volume  data  sets.  It 
may  be  noted  that  the  control  mesh  describing  the  smooth  surfaces  shown  in  Fig.8.8 
has  only  384  triangular  faces  with  a  total  of  194  vertices  as  mentioned  earlier.  This 
experiment  clearly  shows  that  our  model  can  be  used  to  track  a  shape  of  interest 
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Figure  8.8.  Snapshots  from  the  animation  of  canine  heart  motion  over  a  cardiac  cycle 
using  the  dynamic  butterfly  subdivision  model. 
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Figure  8.9.  (a)  The  initialized  model  along  with  the  data  set;  (b)  the  fitted  model 
with  two  subdivisions  on  the  initial  mesh,  along  with  attached  springs  for  editing. 
The  model  after  editing  (c)  at  lower  resolution,  (d)  at  the  same  resolution  of  the 
fitted  model,  and  (e)  at  higher  resolution. 


from  a  set  of  time  dependent  volume  data  sets  in  an  efficient  maimer.  Note  that  no 
other  existing  purely  geometric  subdivision  surface  technique  can  be  used  with  (time 
varying)  continuous  data  sets. 

8.4     Multiresolution  Editing 

The  proposed  dynamic  subdivision  techniques  present  hierarchical  shape  recov- 
ery and  shape  modeling  within  a  common  framework,  where  the  modeler  can  scan 
in  3D  data  points  of  a  prototype  model,  recover  the  underlying  shape  from  the  point 
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set,  and  then  edit  the  recovered  shape.  The  multiresolution  representation  of  the 
evolving  control  mesh  as  presented  in  Chapter  6  enables  the  modeler  to  edit  the 
smooth  limit  surface  at  any  desired  level.  Within  the  proposed  multiresolution  dy- 
namic framework,  the  modeler  does  not  need  to  build  a  model  from  scratch  (unlike 
other  shape  modeling  techniques),  and  there  is  no  need  of  using  computationally  in- 
tensive remeshing  techniques  for  multiresolution  representation  (unlike  other  shape 
recovery  techniques).  The  idea  is  illustrated  in  Fig.8.9,  where  the  shape  of  a  head  is 
first  recovered  from  a  range  data  set  and  then  edited  at  various  levels  using  spring 
forces  from  two  points  in  3D.  The  effect  was  more  global  in  the  lowest  level  edited, 
and  the  effect  became  increasingly  local  as  the  level  of  editing  was  increased. 

8.5     Natural  Terrain  Modeling 

In  this  section,  a  novel  technique  for  synthesizing  realistic  terrain  models  from 
elevation  data  is  developed  by  using  the  dynamic  finite  element  method-based  subdi- 
vision surface  model  for  surface  reconstruction,  and  a  variant  of  the  successive  random 
addition  method  for  fractal  surface  synthesis.  Roughness  is  an  essential  characteristic 
of  the  natural  terrains  and  hence  traditional  surface  reconstruction  methods  using 
smoothness  constraints  do  not  yield  the  desired  results.  Fournier  et  al.  [33]  first  pro- 
posed a  random  displacement  technique  for  synthesizing  fractal  surfaces  which  was 
later  modified  by  Saupe  [84]  in  his  successive  random  addition  scheme  of  generating 
fractals.  Yokoya  et  al.  [107]  improved  these  schemes  by  adding  data  constraints. 
Szeliski  and  Terzopoulos  [94]  proposed  constrained  fractal  surfaces  using  a  Gibbs 
sampler  algorithm  which  was  later  improved  by  Vemuri  et  al.   [102,  103].   Arakawa 
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and  Krotkov  [1]  refined  the  original  Gibbs  sampler  technique  by  redefining  the  tem- 
perature (control)  parameter  to  obtain  better  control  of  roughness  in  the  fitted  (con- 
strained) fractal  surface.  However,  the  (CPU)  execution  times  reported  in  their  work 
are  very  high,  thus  making  their  scheme  unattractive  for  many  applications. 

All  the  techniques  for  natural  terrain  modeling  mentioned  above  usually  needs 
a  grid  of  very  large  size  to  model  realistic  terrains,  especially  with  irregularly  spaced 
data.  Riimelin  [82]  developed  an  interesting  fractal  interpolation  algorithm  which 
can  generate  interpolating  fractal  surfaces  for  regularly  or  irregularly  spaced  data. 
However,  this  scheme  is  computationally  intensive  for  large  problems  as  pointed  out 
in  Vemuri  et  al.  [103]. 

The  subdivision  schemes  produce  a  smooth  surface  in  the  limit,  and  hence  they 
are  not  suitable  for  synthesizing  natural  terrains  which  are  rough.  A  scheme  is  devel- 
oped in  this  dissertation  using  a  variant  of  the  successive  random  addition  method 
[84]  such  that  the  limit  surface  looks  like  "natural"  terrain.  In  the  original  successive 
random  addition  technique,  an  equally  spaced  rectangular  grid  is  refined  by  inter- 
polating the  midpoints  of  each  rectangular  cell  (thereby  dividing  each  rectangular 
cell  into  four  rectangular  cells)  and  then  all  grid  positions  are  perturbed  by  addi- 
tion of  a  Gaussian  noise.  This  process  is  carried  out  recursively  to  obtain  a  fractal 
surface  whose  roughness  is  controlled  by  the  variance  of  the  added  Gaussian  noise 
at  different  refinement  levels.  In  this  dissertation,  the  butterfly  subdivision  scheme 
has  been  modified  where  vertex  positions  at  various  levels  of  subdivision  are  perturbed 
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Figure  8.10.  (a)  Discrete  elevation  data  set  (4096  points),  (b)  fitted  dynamic  butterfly 
subdivision  surface  model  with  841  vertices  (without  noise  addition),  and  (c)  fitted 
dynamic  subdivision  surface  model  with  841  vertices  when  Gaussian  noise  is  added. 


by  addition  of  a  Gaussian  noise  whose  variance  controls  the  roughness  of  the  result- 
ing limit  surface.  This  process  of  Gaussian  noise  addition  is  similar  to  that  of  the 
successive  random  addition  method  mentioned  earlier,  the  difference  being  the  ver- 
tex positions  obtained  using  butterfly  subdivision  rules  are  perturbed  instead  of  grid 
points  obtained  through  midpoint  interpolation. 

Two  natural  terrain  synthesis  results  using  dynamic  modified  butterfly  subdivi- 
sion surface  model  are  presented  in  this  section.  The  initialized  model  is  deformed 
by  applying  spring  forces  on  its  limit  surface  from  the  discrete  data  points.  At  each 
time  step,  every  control  vertex  position  is  perturbed  by  adding  a  random  noise  drawn 
from  a  Gaussian  distribution.  The  variance  of  the  Gaussian  distribution  determines 
the  roughness  of  the  synthesized  surface.  In  both  the  examples,  the  initialized  but- 
terfly subdivision  model  has  an  initial  (control)  mesh  with  98  triangular  faces  and 
68  control  vertices.   The  "natural"  looking  limit  surface  of  the  initialized  model  is 
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Figure  8.11.  Synthesized  natural  terrain  of  different  roughness  using  the  dynamic 
butterfly  subdivision  surface  model  with  841  vertices  from  a  data  set  of  3099  elevation 
values. 


deformed  by  the  forces  exerted  from  the  discrete  elevation  data  points.  When  an 
approximate  fit  is  obtained,  the  model  is  subdivided  to  obtain  a  closer  fit  using  more 
degrees  of  freedom  (control  vertices)  of  the  new  initial  mesh.  The  fitted  surface  has 
1568  triangular  faces  and  841  control  vertices  in  both  the  experiments.  It  may  be 
noted  that  synthesis  of  same  quality  natural  terrains  using  the  existing  techniques 
requires  a  large  number  of  grid  points  (of  the  order  of  105)  [1,  82,  94,  103]  and  hence 
the  proposed  technique  provides  a  more  compact  representation  of  the  synthesized 
terrain.  The  elevation  data  values  are  scaled  to  fit  an  unit  cube  and  the  variance 
of  added  noise  is  10-4  for  the  synthesized  fractal  surface  depicted  in  Fig.8.10.  The 
corresponding  value  of  noise  variance  for  fractal  surfaces  depicted  in  Fig.8. 11(a),  (b) 
and  (c)  are  10~6,  10~5  and  10~4  respectively.  In  the  first  experiment,  4096  elevation 
data  points  are  used  whereas  the  second  data  set  comprised  of  3099  elevation  values. 
The  error  in  fit  is  approximately  1%  in  both  the  examples. 


CHAPTER  9 

CONCLUSIONS  AND  FUTURE  WORK 


In  this  chapter,  the  contributions  of  this  dissertation  are  summarized.  Future 
research  directions  of  the  deformable  subdivision  surface  model  are  also  discussed. 

9.1     Conclusions 

In  this  dissertation,  a  dynamic  framework  has  been  developed  for  subdivision 
surfaces.  This  dynamic  framework  enhances  the  applicability  of  subdivision  surfaces 
in  geometric  modeling  applications  where  the  modeler  can  directly  and  intuitively 
manipulate  the  smooth  limit  surface  using  synthesized  forces.  The  proposed  frame- 
work is  very  useful  for  hierarchical  shape  recovery  from  large  range  and  volume  data 
sets,  as  well  as  for  non-rigid  motion  tracking  from  a  temporal  sequence  of  data  sets. 
Multiresolution  representation  of  the  initial  mesh  controlling  the  smooth  limit  sur- 
face enables  global  and  local  editing  of  the  shape  as  desired  by  the  modeler.  A  novel 
technique  for  synthesizing  natural  terrain  models  from  sparse  elevation  data  using  the 
dynamic  framework  in  conjunction  with  a  variant  of  the  successive  random  addition 
technique  has  also  been  presented  in  this  dissertation. 

Several  theoretical  contributions  has  also  been  made  in  this  dissertation.  Local 
parameterization  techniques  suitable  for  embedding  the  geometric  subdivision  surface 
model  in  a  physics-based  modeling  framework  has  been  developed.    Specific  local 

155 


156 


parameterization  techniques  have  been  fully  worked  out  for  Catmull-Clark,  modified 
butterfly  and  Loop  subdivision  schemes.  Techniques  for  assigning  material  properties 
to  geometric  subdivision  surfaces  have  been  derived,  motion  equation  for  the  dynamic 
model  has  been  formulated,  advantages  of  deformable  models  are  incorporated  into 
conventional  subdivision  schemes,  dynamic  hierarchical  control  has  been  introduced, 
multiresolution  representation  of  the  control  mesh  is  derived  and  an  unified  approach 
for  deriving  subdivision  surface-based  finite  elements  has  been  presented. 

The  proposed  dynamic  framework  has  a  promising  future  in  computer  graph- 
ics, geometric  modeling  and  scientific  visualization.  Furthermore,  the  finite  element 
techniques  proposed  in  this  dissertation  should  be  of  great  interest  to  the  engineering 
design  and  analysis  community  as  well. 

9.2     Future  Directions 

The  proposed  dynamic  framework  enhances  the  applicability  of  the  subdivision 
surfaces  in  various  applications,  but  still  there  are  lots  of  research  need  to  be  done 
to  meet  the  growing  demands  of  modeling  and  other  applications.  In  the  rest  of  this 
section,  several  key  issues  needing  attention  are  discussed. 

9-2.1 Automatic  Change  of  Topology 

In  the  current  framework,  the  evolving  model  can  not  change  its  topology  if 
needed.  Automatic  change  of  topology  for  dynamic  subdivision  surface  models  is 
a  very  challenging  and  important  research  topic.  The  topologically  adaptable  de- 
formable models  developed  by  Malladi  et  al.    [54]  and  Mclnerney  and  Terzopoulos 
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[63]  can  provide  some  insight  about  how  to  deal  with  topology  changing  in  deformable 
subdivision  surface  models. 

9.2.2  Local  Refinement 

In  some  applications,  especially  in  shape  recovery,  it  is  more  meaningful  to 
locally  refine  the  smooth  limit  surface,  i.e.,  introduction  of  new  degrees  of  freedom 
only  in  the  regions  where  more  details  need  to  be  recovered.  Currently,  this  refinement 
is  global  where  new  degrees  of  freedom  is  introduced  throughout  the  smooth  limit 
surface.  Local  refinement  implies  distribution  of  degrees  of  freedom  at  a  number  of 
subdivision  levels  instead  of  one  particular  subdivision  level.  This  is  a  very  challenging 
but  important  research  issue. 

9.2.3  Automatic  Model  Parameter  Selection 

The  model  parameters  like  mass,  damping,  rigidity,  bending  are  selected  by  the 
user  in  the  current  implementation.  The  choice  of  force  constants  and  time  step 
is  also  manual.  Numerical  stability  of  the  deforming  model  is  very  sensitive  to  the 
choice  of  the  parameter  values.  Research  needs  to  be  done  on  how  to  select  these 
parameters  automatically  in  different  applications. 

9.2.4  Constraint  Imposition 

In  some  applications,  there  might  be  a  need  to  impose  geometric  as  well  func- 
tional constraints  on  the  dynamic  subdivision  surface  models.  Position  and  normal 
constraints  are  examples  of  possible  constraints.  It  is  somewhat  easier  to  impose  these 
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constraints  on  interpolators  dynamic  subdivision  surface  models,  but  considerable  re- 
search efforts  need  to  be  put  for  constraint  imposition  in  the  case  of  approximating 
dynamic  subdivision  surface  models. 

9.2.5  Recovery  of  Sharp  Features 

Recovery  of  sharp  features  using  dynamic  subdivision  surface  models  is  an  open 
research  problem.  Even  though  piecewise  smooth  subdivision  surfaces  [40]  has  been 
developed,  incorporating  sharp  features  in  a  physically  meaningful  way  for  an  evolving 
subdivision  surface  model  is  a  challenging  research  issue. 

9.2.6  Automatic  Model  Initialization 

Currently  model  initialization  is  interactive.  However,  automatic  model  initial- 
ization techniques  have  been  recently  proposed  by  Delingette  [24].  Similar  techniques 
for  deformable  subdivision  surface  initialization  need  to  be  developed. 

9.2.7  Improved  Synthesized  Force  Fields 

Recently  Xu  and  Prince  [106]  have  proposed  new  type  offeree  fields  which  makes 
the  model  insensitive  to  the  initialization.  It  would  be  nice  to  incorporate  these  force 
fields  in  the  current  system  for  better  performance  in  the  model  fitting  context. 
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